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CHAPTER  I 


SUMMARY 


This  contract  was  used  to  support  work  on  several  different  problems.  Each  chapter,  below, 
describes  a  different  problem.  In  this  chapter,  we  summarize  the  results  of  those  later  chapters. 


SUMMARY  OF  CHAPTER  II 


We  examined  the  theoretical  effects  of  air  pressure  loading  on  surface  displacements  and  surface 
gravity.  We  concluded  that  global  baseline  changes  due  to  this  effect  can  frequently  be  2  cm  or  more, 
and  that  gravity  perturbations  can  be  as  large  as  3-6  microgals.  The  effects  at  a  station  can  be  ade¬ 
quately  calculated  using  real  pressure  data  within  about  1000  km  of  that  station.  For  inland  points,  ade¬ 
quate  results  can  be  obtained  using  only  two  parameters  to  describe  the  pressure,  as  described  by  Rab- 
bel  and  Zschau  (J.Geophys.,  56,  81-99,  1985).  At  points  near  the  coast,  however,  more  complete  pres¬ 
sure  data  are  required. 


SUMMARY  OF  CHAPTER  III 


We  investigated  the  effects  of  core  structure  on  the  forced  nutations.  VLBI  observations  (from 
Herring,  et  ai„  and  Gwinn,  ct  al.)  and  tidal  gravity  studies  (from  Levine,  ct  al.,  and  Neubcrg,  ct  al.) 
suggest  that  the  fluid  core  resonance  in  the  nutations  and  tides  should  have  an  inertial  space  period  of 
around  433  days,  instead  of  the  460  days  predicted  for  a  hydrostatically  pre-stressed  earth.  We 
developed  a  theoretical  model  of  the  resonance  for  an  earth  which  is  non-hydrostatically  pre-stressed, 
but  which  has  a  fluid  outer  core.  We  find,  for  this  model,  that  the  resonance  is  sensitive  to  the  1=2, 
m=0  spherical  harmonic  term  in  the  non-hydrostatic  structure,  but  not  to  any  other  spherical  harmonic 
component. 


SUMMARY  OF  CHAPTER  IV 


We  found  that  there  can  be  lateral,  non-hydrostatic  structure  inside  th.  id  core,  caused  by  grav¬ 
itational  forcing  from  the  mantle,  from  the  inner  core,  or  from  topography  on  the  corc/mantlc  or  inner 
core/outer  core  boundaries.  The  possibility  of  such  structure  has  been  overlooked  in  the  past,  and  it  can 
have  important  implications  for  geodetic  and  seismic  studies  of  the  core.  We  developed  the  differential 
equations  which  describe  the  structure,  and  considered  the  effects  of  the  solutions  on  nutation  and  tidal 
constraints  for  core/mantlc  boundary  topography.  We  concluded  that  those  constraints  could  accommo¬ 
date  large  1=2,  m=0  topography  on  the  boundary,  if  there  was  a  thin,  low  density  fluid  layer  just 
bui.caih  !h°  boundary. 
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SUMMARY  OF  CHAPTER  V 


We  began  addressing  the  problem  of  the  ocean’s  response  to  pressure.  It  is  important  to  under¬ 
stand  this  response  when:  (a)  modelling  the  pressure-induced  and  wobble-induced  crustal  deformation  at 
gra'd'y  VI. BI,  satellite  ranging,  etc.,  sites  within  a  few  hundred  km  of  the  coast;  (b)  estimating  the 
effects  of  global  pressure  on  polar  motion;  (c)  using  satellite  altimeter  data  to  study  the  wind-driven  cir¬ 
culation  of  the  ocean.  We  analyzed  GEOSAT  altimeter  data  to  study  this  response.  (GEOSAT  data  are 
particularly  helpful  to  investigate  the  short  period  response;  tide  gauge  data  are  helpful  to  study  the 
long  period  response  -  see  item  (6),  below.) 

Our  first  work  in  this  area  has  been  to  assess  the  reliability  of  global  pressure  data.  We  found 
substantial  differences  between  results  from  different  meteorological  centers,  particularly  over  the 
Southern  Pacific. 


SUMMARY  OF  CHAPTER  VI 

We  used  monthly  tide  gauge  data  from  the  Permanent  Service  for  Mean  Sea  Level  in  Bidston, 
England,  to  investigate  certain  long-period  disturbances  in  the  ocean.  We  found;  (a)  the  18.6  year  lunar 
tide  and  14  month  pole  tide  are  close  to  equilibrium  in  the  global  ocean  (it  is  important  to  understand 
these  ocean  tides  when  interpreting:  (i)satcllitc  ranging  results  for  the  18.6  year  earth  tide,  and  (ii)  polar 
motion  observations  of  the  Chandler  wobble  period,  in  terms  of  mantle  anclasticity);  (b)  the  response  of 
the  ocean  to  pressure  is  very  close  to  inverted  barometer  at  periods  of  two  months  and  longer;  (c)  the 
global  rise  in  sea  level  over  the  last  80  years,  after  removing  the  effects  of  post-glacial  rebound  on  indi¬ 
vidual  tide  gauge  records,  has  been  between  1.1  mm/yr  and  1.9  mm/yr. 


CHAPTER  II 


DISPLACEMENTS  OF  THE  EARTH’S  SURFACE  DUE  TO  ATMOSPHERIC  LOADING: 
EFFECTS  ON  GRAVITY  AND  BASELINE  MEASUREMENTS 


SUMMARY 


Atmospheric  mass  loads  and  deforms  the  earth’s  crust.  By  performing  a  convolution  sum 
between  daily,  global  barometric  pressure  data  and  mass  loading  Green’s  functions,  we  estimate  the 
time  dependent  effects  of  atmospheric  loading,  including  those  associated  with  short-term  synoptic 
storms,  on  surface  point  positioning  measurements  and  surface  gravity  observations.  We  calculate  the 
response  for  both  an  oceanless  earth  and  an  earth  with  an  inverted  barometer  ocean.  Load  responses  for 
near-coastal  stations  are  significantly  affected  by  the  inclusion  of  an  inverted  barometer  ocean.  Peak  to 
peak  vertical  displacements  are  frequently  15-20  mm  with  accompanying  gravity  perturbations  of  3-6 
pGal.  Baseline  changes  can  be  as  large  as  20  mm  or  more.  The  perturbations  are  largest  at  higher  lati¬ 
tudes  and  during  winter  months.  These  amplitudes  are  consistent  with  the  results  of  Rabbel  and  Zschau 
[1985],  who  modeled  synoptic  pressure  disturbances  as  Gaussian  functions  of  radius  around  a  central 
point.  Deformation  can  be  adequately  computed  using  real  pressure  data  from  points  within  about  1000 
km  of  the  station.  Knowledge  of  local  pressure,  alone,  is  not  sufficient.  Rabbel  and  Zschau’s 
hypothesized  corrections  for  these  displacements,  which  use  local  pressure  and  the  regionally  averaged 
pressure,  prove  accurate  at  points  well  inland  but  are,  in  general,  inadequate  within  a  few  hundred 
kilometers  of  the  coast. 


INTRODUCTION 

Recent  improvements  in  geodetic  instrumentation  and  measurement  techniques  have  made  the 
detection  of  crustal  movements  over  short  time  periods  feasible.  It  is  expected  that  space  techniques 
such  as  very  long  baseline  interferometry  (VLBI)  and  satellite  laser  ranging  will  soon  have  the  capabil¬ 
ity  of  measuring  distances  over  long  baselines  to  within  a  centimeter  or  better  (see,  for  example,  Her¬ 
ring  [1986],  and  Cohen  and  Smith  [1985].  Relative  gravity  observations  used  to  detect  vertical  crustal 
motions  have  already  achieved  sub-microGal  accuracies  over  long  time  periods  (see  for  example, 
Richter  [1983],  Levine  et  al.  [1986],  and  Goodkind  [1986]).  These  techniques  promise  to  provide  valu¬ 
able  information  about  tectonic  deformation  of  the  earth.  Since  nontectonic  processes  can  also  cause 
crustal  motions,  their  effects  need  to  be  considered  as  well  and,  if  necessary,  removed  from  the  data.  In 
this  chapter  we  consider  displacements  caused  by  global  barometric  pressure  fluctuations  around  the 
earth. 

Variations  in  the  horizontal  distribution  of  atmospheric  mass  induce  deformations  within  the 
earth.  The  interaction  between  the  atmosphere  and  the  earth  is  through  pressure  loading  at  the  surface 
and,  particularly  at  longer  wavelengths,  through  the  gravitational  attraction  of  the  atmospheric  mass. 
The  displacement  of  the  surface  is  primarily  vertical,  and  there  is  a  change  in  the  acceleration  of  grav¬ 
ity  at  the  displaced  surface. 

Horizontal  changes  in  atmospheric  mass  can  be  inferred  from  surface  atmospheric  pressure  meas¬ 
urements.  Stolz  and  Larden  [1979]  concluded  that  global  seasonal  fluctuations  in  barometric  pressure 
would,  in  general,  contribute  less  than  a  centimeter  to  surface  displacements.  The  largest  pressure  vari¬ 
ations,  however,  are  those  associated  with  synoptic  scale  storms.  Because  the  average  life  span  of  these 
storms  is  approximately  5  days  [Klein,  1956],  they  do  not  contribute  appreciably  to  the  seasonal  varia¬ 
bility  in  global  pressure. 
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A  simple  calculation  using  elastic  Green’s  functions  for  the  spherically  symmetric  Gutenberg- 
Bullen  A  earth  model  from  Farrell  [1972]  shows  that  a  uniform  pressure  increase  of  20  mbar  over  a 
disk  of  diameter  2000  km  causes  a  depression  of  10  mm  at  the  center  of  the  disk.  This  displacement 
scales  linearly  with  the  pressure  over  the  disk  and  also  increases,  although  more  slowly  than  linearly, 
with  the  diameter  of  the  disk.  For  example,  if  the  20-mbar  pressure  is  distributed  over  a  disc  with  only 
a  1000-km  diameter,  the  central  depression  is  about  6  mm.  Since  large  synoptic  storms  can  have  spa¬ 
tial  scales  of  about  2000  km  and  peak-to-peak  pressure  fluctuations  which  may  be  as  large  as  40-50 
mbar,  it  is  probable  that  vertical  displacements  could  frequently  be  10  mm  or  larger  along  a  storm 
track. 

Rabbel  and  Zschau  [1985]  estimated  the  deformation  caused  by  these  synoptic  storms  by  assum¬ 
ing  that  the  accompanying  pressure  disturbance  is  a  Gaussian  function  of  radius  about  some  central  high 
or  low  and  then  convolving  that  pressure  distribution  with  Farrell’s  Green’s  functions.  Their  results 
suggest  that  vertical  displacements  of  the  earth’s  surface  of  10-20  mm  can  occur.  They  found  that  they 
could  adequately  predict  the  deformation  caused  by  a  Gaussian  pressure  distribution  by  using  a  simple 
regression  equation  together  with  estimates  of  the  local  pressure  at  the  point  of  interest  and  of  the  aver¬ 
age  surface  pressure  with  a  2000-km  radius  of  the  point. 

In  this  chapter  we  extend  and  assess  Rabbel  and  Zschau’s  results  by  convolving  Farrell’s  Green’s 
functions  with  real.  Interpolated  global  pressure  data.  Radial  and  horizontal  displacements  as  well  as 
changes  in  the  gravitational  acceleration  for  various  stations  are  calculated.  Baseline  changes  between 
stations  are  also  examined. 

In  particular,  we  address  the  following  questions.  First,  are  the  displacements  predicted  using  real 
data  as  potentially  important  as  suggested  by  the  results  of  Rabbel  and  Zschau?  We  conclude  that  they 
are.  We  find  that  displacements  are  generally  larger  at  higher  latitudes,  and  during  the  winter  months. 
Typically,  peak-to-peak  displacements  of  a  centimeter  or  two  occur  during  each  winter  month  at  higher 
latitudes. 

Second,  a  change  in  atmospheric  pressure  can  cause  a  change  in  sea  level  height  which  further 
deforms  the  earth.  How  important  is  it  to  include  this  ocean-induced  deformation  in  the  calculations? 
We  find  that  if  the  ocean  responds  to  pressure  variations  as  an  inverted  barometer,  the  effect  on  vertical 
displacements  can  be  greater  than  10  mm  at  stations  on  the  coast,  and  often  10  mm  within  several  hun¬ 
dred  kilometers  of  the  coast. 

Third,  do  Rabbel  and  Zschau's  regression  equations  provide  adequate  approximations  to  the 
deformation  predicted  from  real  data?  We  find  that  the  equations  work  well  at  inland  points  but  are,  in 
general,  inadequate  at  points  within  a  few  hundred  kilometers  of  the  coast,  assuming  that  an  inverted 
barometer  ocean  response  is  appropriate. 


METHOD 

Farrell  [1972]  constructed  Green’s  functions  to  describe  the  response  of  an  elastic  earth  to  a  point 
load  on  its  surface.  Gravitational  accelerations,  displacements,  lilts,  and  strains  can  be  determined  for 
any  surface  load  by  evaluating  a  convolution  sum  between  the  Green’s  functions  and  the  load  distribu¬ 
tion.  The  effects  of  anelasticity  at  the  periods  of  interest  here  are  almost  certain  to  be  less  than  \d% 
and  so,  judging  by  the  results  below,  can  be  ignored. 

We  obtained  (from  the  National  Center  for  Atmospheric  Research  in  Boulder,  Colorado)  the 
National  Meteorological  Center  (NMC)  global  atmospheric  pressure  values  from  1979  through  1984. 
The  NMC  calculated  surface  pressure  values  twice  daily  for  a  2.5°  x  2.5°  grid  over  the  entire  surface  of 
the  earth.  Because,  we  are  interested  only  in  time  dependent  displacements,  we  subtract  the  atmos¬ 
pheric  pressure  for  an  arbitrary  day  from  all  other  daily  values. 

To  estimate  the  displacement  at  a  surface  point,  we  partition  the  earth  into  a  similar  2.5°  x  2.5° 
grid  system.  Grid  elements  that  lie  within  5°  of  the  field  point  are  further  subdivided  into  smaller 
units.  Values  for  the  Green’s  functions  are  assigned  to  each  grid  element  and  arc  determined  by  the 
angular  distance  between  the  element  and  the  surface  point.  The  Green’s  functions  arc  multiplied  by 
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the  daily  atmospheric  load,  and  the  results  are  summed.  In  this  way,  we  can  determine  the  radial  and 
horizontal  displacements  and  the  change  in  the  gravitational  acceleration  at  any  point  on  the  surface  of 
the  earth.  We  calculate  the  change  in  gravity  due  to  the  elastic  deformation  of  the  solid  earth  plus  that 
due  to  the  gravitational  attraction  of  all  displaced  atmospheric  mass  except  for  the  excess  mass  directly 
above  the  point.  The  contribution  from  that  excess  mass  can  be  crudely  approximated  by  multiplying 
the  local  pressure  in  millibars  by  the  factor  0.42  pGal/mbar,  determined  by  assuming  the  excess  mass  is 
an  infinite  plane  of  uniform  density.  This  latter  pressure  contribution  is  typically  larger  than  the  effects 
we  are  considering  here  by  a  factor  of  3  and  is  already  routinely  included  in  many  data  analyses.  The 
effects  of  the  mass  directly  overhead  can  not  be  exactly  described  with  a  constant  factor  of 
0.42  pGal/mbar,  since  the  mass  is  not  a  uniform  plane.  A  better  way  to  remove  those  effects  is  to  fit 
them  to  the  data  [Warburton  and  Goodkind ,  1977;  Spratt,  1982;  Levine  et  at.,  1986]. 


RESULTS 

Deformations  and  gravitational  accelerations  are  calculated  using  both  an  oceanless  earth  model 
and  a  model  in  which  the  oceans  respond  as  an  inverted  barometer  to  atmospheric  loading.  The 
inverted  barometer  hypothesis  states  that  for  every  millibar  increase  in  atmospheric  pressure,  the  sea 
surface  depresses  locally  by  1  cm.  This  implies  that  the  ocean  basins  experience  none  of  the  forces 
associated  with  barometric  pressure  fluctuations  and  is  equivalent  to  setting  the  total  incremental  mass 
load  over  the  ocean  basins  equal  to  zero.  However,  to  conserve  oceanic  mass,  the  inverted  barometer 
hypothesis  must  be  modified  slightly.  If  there  is  a  net  increase  or  decrease  in  the  mass  of  air  above  the 
oceans,  the  seafloor  experiences  a  uniform  pressure  D  acting  everywhere  on  the  earth’s  surface  under 
the  oceans  and  equal  to 
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where  AP  is  the  local  change  in  pressure,  A  is  the  surface  area  of  the  oceans,  and  the  integral  is  over 
the  ocean  surface.  Results  indicate  that  the  effects  of  the  pressure  D  on  vertical  displacements  are  no 
more  than  3  mm  for  coastal  geodetic  stations  and  are  even  smaller  at  inland  points. 

It  is  likely  that  the  ocean  does  respond  to  pressure  as  an  inverted  barometer  at  periods  of  weeks 
to  a  few  years.  Wunsch  [1972]  compared  local  pressure  and  tide  gauge  data  taken  at  Bermuda  and  con¬ 
cluded  that  a  local  inverted  barometer  response  was  probably  appropriate  at  periods  as  short  as  3  or  4 
days.  On  the  other  hand,  the  decidedly  nonequilibrium  diurnal  ocean  tides  imply  that  the  global 
response  is  certainly  not  an  inverted  barometer  at  periods  close  to  1  day.  Furthermore,  Eubanks  et  al. 
[1985]  have  concluded  from  recent  polar  motion  data  that  the  southern  ocean  may  differ  significantly 
from  the  inverted  barometer  response.  We  are  thus  not  certain  that  our  global  inverted  barometer 
assumption  is  appropriate  for  the  short-period  synoptic  storms  which  most  concern  us.  (See  Chelton 
and  Enfield  [1986]  for  a  summary  of  observational  evidence  for  the  inverted  barometer  response.) 

We  attempt  to  bound  the  actual  atmospheric-induced  deformations  by  using  the  oceanless  earth 
and  the  inverted  barometer  model  as  end  points.  The  larger  estimates  of  atmospheric  loading  effects 
are  obtained  for  the  oceanless  earth  model.  Figure  II.  1  is  a  graph  of  the  calculated  surface  displace¬ 
ments  and  associated  gravity  perturbations  at  Boulder,  Colorado,  during  1980  for  such  a  model.  The 
largest  peak-to-peak  barometric  pressure  fluctuations  (Figure  Il.la)  and  crustal  responses  occur  during 
the  winter  months  where  radial  displacements  often  exceed  10  mm  (Figure  11.16).  Horizontal  displace¬ 
ments  are  also  calculated  (not  shown)  but  are  found  to  be  only  between  1  and  2  mm.  Fluctuations  in 
the  local  value  of  gravity  (besides  those  approximately  described  by  the  0.42  pGal/mbar  factor  dis¬ 
cussed  above)  are  shown  in  Figure  II. lc.  Perturbations  are  determined  to  be  about  2-4  pGal  and  are 
mostly  due  to  the  vertical  displacement  of  the  surface.  The  calculated  displacements  and  gravity 
fluctuations  are  somewhat  smaller  during  the  summer  months,  with  maximum  peak  to  peak  redial  dis¬ 
placements  of  about  10  mm. 
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At  higher  latitudes  the  loading  effects  are  systematically  larger  due  to  the  larger  pressure  varia¬ 
tions  there.  For  example,  at  Onsala,  Sweden,  during  1980  our  calculations  for  an  oceanless  earth 
predict  that  vertical  crustal  motions  of  20-30  mm  (Figure  11.26)  and  gravity  perturba'ions  of  6-8  pGal 
(Figure  U.2„)  occur  regularly  during  the  winter  months.  Horizontal  displacements  are  larger  as  well  but 
are  still  small,  having  maximum  magnitudes  of  2-3  mm. 

We  then  recompute  the  load  responses  using  the  inverted  barometer  assumption.  The  addition  of 
oceans  to  the  earth  model  has  a  greater  effect  on  the  points  closer  to  the  coast.  Radial  displacements  at 
Onsala,  situated  on  the  North  Sea  and  close  to  the  North  Atlantic  Ocean,  are  reduced  by  as  much  as  10 
mm  (Figure  11.36)  and  gravity  perturbations  are  decreased  by  up  to  2  pGal.  A  large  reduction  is  not 
surprising  since  when  the  inverted  barometer  assumption  is  applied,  the  crust  under  the  nearby  ocean  no 
longer  experiences  the  effects  of  large  pressure  fluctuations  and  therefore  no  longer  contributes  to  the 
deformations  at  Onsala.  The  load  response  at  inland  stations  is  less  affected  by  the  addition  of  an 
inverted  barometer  ocean  to  the  model.  For  example,  at  Boulder  the  disp'acemcnts  calculated  for  an 
inverted  barometer  model  differ  by  only  a  couple  of  millimeters  from  those  for  the  oceanless  earth 
model  (Figure  II. 3a).  Still,  we  find  that  differences  of  at  least  a  centimeter  can  persist  up  to  several 
kilometers  inland. 

The  results  of  our  calculated  perturbations  in  gravity  were  used  by  Levine  el  al.  [1986]  to  correct 
gravity  measurements  for  atmospheric  effects.  They  used  a  LaCoste  and  Romberg  gravimeter  to  record 
continuously  the  local  changes  in  the  gravitational  acceleration  at  Boulder  from  May  through  December 
1983.  By  removing  the  gravitational  effects  due  to  global  pressure  variations,  Levine  et  al.  were  able  to 
reduce  the  amplitude  of  their  residuals  by  about  one  third  and  improve  the  agreement  between  calcu¬ 
lated  and  observed  tidal  admittances. 

Baseline  lengths  are  affected  by  displacements  of  the  end  points.  Because  the  displacements  are 
mostly  radial,  the  effects  are  largest  on  very  long  baselines.  As  an  example,  the  effects  on  the  Onsala 
to  Haystack  Observatory  baseline  during  1984  are  shown  in  Figure  II.4a.  This  baseline  is  too  short 
(only  about  one  earth  radius)  to  maximize  the  effect,  but  it  is  one  of  the  more  carefully  and  frequently 
measured  VLBI  baselines  [Herring  et  ai.  1986], 

Vertical  motions  at  both  Onsala  and  Haystack  can  be  as  large  as  30  mm  for  the  oceanlcss  earth 
model.  Corresponding  baseline  changes  of  15-20  mm  arc  common,  and  maximums  of  35  mm  can 
occur  during  the  winter.  Baseline  changes  are  decreased  by  using  the  inverted  barometer  assumption 
(Figure  11.46),  where  maximum  variations  are  of  the  order  of  20  mm.  We  find  similar  baseline  changes 
for  other  medium  to  long  baselines  at  high  latitudes. 

The  effects  on  a  baseline  depend  on  the  relative  motion  of  the  end  points.  An  examination  of  the 
vertical  motions  at  Onsala  and  Haystack  indicates  that  the  earth  deforms  approximately  in-phase  at 
these  points,  thus  increasing  the  atmospheric  effects  on  baseline  changes.  When  pressure  variations 
associated  with  synoptic  storms  around  the  globe  are  expanded  in  terms  of  spherical  harmonics,  T". 
the  largest  coefficients  are  found  to  occur  at  n= 6  and  7  [Sprait,  1982],  These  harmonics  correspond  to 
a  wavelength  of  around  6000  km,  which  is  about  equal  to  the  distance  over  the  earth’s  surface  between 
Onsala  and  Haystack,  a  result  consistent  with  the  in-phase  deformation  at  the  two  stations. 

Herring  et  al.  [1986]  measured  the  Onsala-Haystack  baseline  with  VLBI  on  34  different  days 
during  1980-1984.  We  compared  our  estimates  for  this  baseline  with  their  results  and  found  no 
significant  correlation.  However,  our  estimates  with  the  inverted  barometer  ocean  for  those  particular 
days  arc  mostly  within  only  a  few  millimeters  of  each  other,  with  only  one  outlying  point  as  large  as  10 
mm.  By  contrast,  errors  in  'he  VLBI  results  due  to  uncertainties  in  the  atmospheric  water  vapor  are 
probably  around  20  mm  or  so  for  this  data  (T.  Herring,  personal  communication,  1986). 

The  various  results  described  above  are  consistent  in  amplitude  with  the  general  results  of  Rabbel 
and  Zschau  [1985].  They  reinforce  Rabbel  and  Zschau’s  conclusion  that  atmospheric  pressure  correc¬ 
tions  to  geodetic  point  positioning  and  gravity  measurements  should  be  considered  in  future  analyses. 

We  next  estimate  the  spatial  extent  of  real  data  required  in  order  to  accurately  apply  atmospheric 
corrections  using  the  method  described  in  this  chapter.  Dotted  curve  in  Figure  11.5  shows  the  predicted 
radial  displacements  at  Boulder  obtained  by  considering  only  those  pressure  changes  in  the  NMC  data 
which  occur  at  points  within  a  10°  (1000  km  radius)  of  Boulder.  Solid  curve  shows  the  results  obtained 
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using  the  entire  global  data  set.  The  difference  between  these  results  and  those  obtained  using  all  the 
data  is  always  less  than  3  mm,  suggesting  that  it  is  not  presently  necessary  to  consider  the  entire  atmos¬ 
pheric  load.  However,  we  also  find  that  it  is  not  sufficient  to  use  only  the  pressure  variations  close  to 
the  point  of  interest.  For  example,  calculations  which  include  only  the  pressure  fluctuations  within  1° 
(about  100  km)  of  Boulder  can  differ  from  the  global  results  by  as  much  as  10  mm  (dashed  curve  in 
Figure  11.5). 

To  explore  this  result  further,  we  calculate  the  vertical  displacements  at  Boulder  on  January  8, 
1980,  using  several  different  cutoff  angles  for  the  data.  The  vertical  displacement  computed  using  the 
entire  global  data  set  is  large  on  this  day,  roughly  13  mm.  The  displacements  predicted  using  different 
cutoff  angles  are  plotted  in  Figure  11.6.  The  differences  from  the  global  result  are  5  mm  for  a  5°  cutoff, 
2.5  mm  for  a  10°  cutoff,  and  less  than  0.5  mm  for  a  20°  cutoff.  Differences  can  be  somewhat  larger  on 
other  days  and  at  higher  latitudes,  since  the  maximum  displacements  are  larger  there.  These  results 
suggest  that  pressure  data  from  at  least  the  surrounding  1000  km  is  desirable,  given  the  accuracy  levels 
of  the  various  new  geodetic  techniques.  This  distance  is  consistent  with  the  observation  that  the  largest 
storms  are  on  the  order  of  2000  km  across. 

It  would  be  impractical  if  a  barometer  array  had  to  be  installed  whenever  an  observation  was 
made,  to  determine  pressure  variations  within  the  surrounding  1000  km.  Fortunately,  interpolated  glo¬ 
bal  pressure  data  are  available  from  the  NMC  with  only  a  few  weeks  delay.  Still,  it  might  be  con¬ 
venient  in  some  cases  if  a  linear  relationship  between  the  local  and  regionally  averaged  pressures  and 
the  local  crustal  displacement  could  be  found.  By  mathematically  modeling  the  pressure  distribution  of 
synoptic  scale  storms  as  an  exponentially  decaying  function  with  radius  from  the  center  of  the  storm 
and  convolving  the  result  with  Farrell’s  Green's  functions.  Rabbet  and  Zschau  [1985]  hypothesized  a 
relationship  between  the  local  pressure,  the  average  pressure  within  2000  km  of  the  station,  and  the 
local  radial  displacement  and  change  in  the  acceleration  of  gravity: 

Radial  displacement  [Rabbet  and  Zschau,  1985] 


U  =  -0.90F,  -  0.35(P,  -  P,) 


Gravity,  besides  the  0.42  pGal/mbar  factor  described  above  [ Rabbet  and  Zschau,  1985], 

G  =  -0.17/*,  -  0.08(P,  -  P,) 

where  U  is  in  millimeters,  G  is  in  microGal,  and  the  pressure  values  Pt  and  P,  are  in  millibars.  Here, 
Ps  is  the  local  pressure  change,  and  Pt  is  obtained  by  averaging  the  pressure  variations  within  a  sur¬ 
rounding  radius  of  2000  km.  Rabbel  and  Zschau  suggested  that  for  an  inverted  barometer  ocean,  P, 
should  be  calculated  by  setting  the  pressure  changes  for  the  ocean  areas  to  zero.  To  use  these  correc¬ 
tions,  one  would  only  need  to  know  the  local  pressure  changes  and  to  use  a  weather  map  to  obtain  a 
value  for  the  average  pressure  in  the  surrounding  area. 

To  assess  the  accuracy  of  these  equations,  we  use  them  to  estimate  vertical  displacements  at  vari¬ 
ous  points,  using  both  the  oceanless  earth  and  the  inverted  barometer  ocean  models.  The  results  are 
compared  with  displacements  calculated  using  the  real  NMC  pressure  data  (for  comparison,  we  let  D 
equal  zero  when  using  the  inverted  barometer  model).  For  an  oceanless  earth  model  and  at  points  well 
inland  for  the  inverted  barometer  model,  the  results  generally  agree  to  within  90%,  suggesting  that  the 
equations  of  Rabbel  and  Zschau  are  useful  for  calculating  atmospheric  corrections  for  inland  points. 
This  is  mostly  confirmation  that  the  pressure  disturbances  associated  with  large  storms  arc,  indeed,  rea¬ 
sonably  well  represented  by  Gaussian  functions.  Contour  plots  of  the  pressure  data  during  times  of 
maximum  displacement  further  reinforce  this  assumption. 

However,  the  agreement  is  not  as  encouraging  near  coasts  for  the  inverted  barometer  ocean.  For 
example,  the  two  sets  of  results  for  Onsala  during  the  winter  of  1980  have  peak  to  peak  discrepancies 
of  as  much  as  10  mm  (Figure  11.7).  Discrepancies  as  large  as  5  mm  or  more  can  persist  for  several 
hundred  kilometers  inland.  The  reason  is  that  when  the  pressure  is  set  equal  to  zero  over  the  oceans, 
the  Gaussian  character  of  the  storms  is  destroyed  sufficiently  that  Rabbel  and  Zschau’s  regression 
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equations  are  no  longer  adequate.  W.  Rabbel  (personal  communication,  1986)  suggests  that  corrections 
for  near-coastal  stations  might  conceivably  be  obtained  by  extending  the  regression  equations  to  include 
a  modification  factor  appropriate  to  each  station.  On  the  other  hand,  since  there  is  only  a  few  weeks 
delay  in  obtaining  the  NMC  sea  level  pressure  data,  accurate  atmospheric  load  corrections  can  be 
obtained  relatively  quickly  for  both  coastal  and  inland  stations  by  performing  a  convolution  sum  with 
the  data  as  outlined  in  this  chapter. 


CONCLUSIONS 

The  magnitude  of  pressure-induced  deformations  of  ihe  earth’s  surface  and  of  the  associated  grav¬ 
ity  changes  are  within  the  accuracies  obtainable  by  present  measurement  techniques.  Using  real  atmos¬ 
pheric  data,  we  conclude  that  vertical  surface  displacements  of  20  mm  and  gravity  perturbations  of 
4  pGal  are  common.  The  displacements  are  generally  the  largest  during  winter  and  at  higher  latitudes. 
These  amplitudes  are  consistent  with  the  conclusions  of  Rabbel  and  Zschau  [1985],  who  used  an  analyt¬ 
ical,  Gaussian  model  for  storm-related  pressure  disturbances.  Rabbel  and  Zschau  proposed  simple 
two-parameter  equations  'o  approximate  these  deformations.  We  find  that  the  equations  are  accurate 
for  points  well  inland  but  near  the  coasts  they  can  differ  from  results  computed  with  real  data  by  as 
much  as  50%.  The  gravity  amplitudes  are  also  consistent  with  the  results  estimated  by  Spratt  [1982] 
using  atmospheric  pressure  data.  (Spratt  used  a  somewhat  different  computational  approach  than  s.) 

At  coastal  stations  and  at  points  within  a  few  hundred  kilometers  of  the  coast,  the  value  of  the 
response  depends  on  the  ocean  model  applied.  At  Onsala,  Sweden,  during  1980,  for  example,  the  max¬ 
imum  vertical  displacements  are  reduced  from  45  mm  for  an  oceanless  earth  to  30  mm  for  an  inverted 
barometer  ocean.  Calculated  maximum  changes  for  the  Onsala- Haystack  baseline  arc  between  35  and 
20  mm  depending  on  the  ocean  model.  Observation  points  more  than  a  few  hundred  kilometers  inland 
are  not  appreciably  affected  by  the  choice  of  ocean  model.  It  is  not  certain  that  the  inverted  barometer 
assumption  is  justified  at  the  time  scale  of  synoptic  storms,  but  we  have  attempted  to  bound  the  true 
surface  displacements  between  the  oceanless  earth  and  inverted  barometer  assumptions.  We  are 
presentlv  working  to  understand  the  global  oceanic  response  to  pressure  in  more  detail. 

The  results  presented  here  indicate  that  air  pressure  corrections  are  large  enough  to  affect 
significantly  many  point  positioning  measurements.  A,Uq"ate  corrections  are  not  possible  if  only  local 
pressure  fluctuations  are  included.  Wc  have  determined  empirically  that  pressure  data  from  a  1000-km 
radius  around  the  observation  point  are  required  to  match  accurately  the  displacement  calculated  from 
global  pressure  data. 
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CHAPTER  II  FIGURE  CAPTIONS 


Figure  II. 1.  (a)  Local  pressure  variations  at  Boulder,  Colorado,  for  winter  1980.  (£>)  Radial  displacements, 
(c)  Gravity  perturbations  (does  not  include  the  0.42  |J.GaI/mbar  factor  discussed  in  text)  Note  that  the  zero 
lines  in  these  plots  have  has  no  physical  significance  since  daily  pressure  was  differenced  with  that  for  an 
arbitrary  day. 

Figure  11,2.  Same  as  Figure  1  for  Onsala,  Sweden,  for  winter  1980. 

Figure  II.3.  (a)  The  solid  line  represents  the  radial  displacements  at  Boulder  for  the  inverted  barometer 
model;  the  dashed  line  represents  the  radial  displacements  for  an  oceanless  earth  model,  (b)  Same  as  Fig¬ 
ure  3 a  but  for  Onsala. 

Figure  II.4.  (a)  Baseline  change  between  Haystack,  Massachusetts  and  Onsala  for  the  oceanless  earth 
model  [  1984],  ( b )  Same  as  Figure  4 a  for  the  inverted  barometer  model. 

Figure  II.5.  Radial  displacements  at  Boulder  for  January  1982  with  the  inverted  barometer  model.  Solid 
curve,  using  data  from  the  entire  earth;  dotted  curve,  considering  only  pressure  changes  within  1000  km  of 
Boulder;  dashed  curve,  considering  only  pressure  changes  within  100  km  of  Boulder. 

Figure  II.6.  Radial  displacements  at  Boulder,  computed  for  January  8, 1980,  using  different  cutoff  angles  0 
(pressure  values  are  set  to  zero  at  all  grid  points  further  than  0  away  from  Boulder),  the  vertical  subsidence 
at  Boulder  computed  using  the  entire,  global  data  set  is  13.4  mm. 

Figure  II.7.  Radial  displacements  at  Onsala  during  winter  1980,  including  inverted  barometer  oceans.  The 
solid  curve  represents  displacements  calculated  from  real  data,  as  discussed  in  this  chapter.  The  dashed 
curve  represents  displacements  calculated  from  Rabbel  and  Zschau’s  regression  equation.  Agreement  is 
much  better  at  stations  a  few  hundred  kilometers  inland. 
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CHAPTER  III 


THE  EARTH’S  FORCED  NUTATIONS:  GEOPHYSICAL  IMPLICATIONS 


INTRODUCTION 

The  Earth’s  nutational  motion  consists  of  periodic  tipping  of  the  Earth  in  space  and  is  caused  by 
the  gravitational  attraction  of  the  Sun  and  Moon.  The  motion  occurs  at  discrete  frequencies  of,  as  seen 
from  the  Earth,  one  cycle  per  day  modulated  by  the  orbital  frequencies  of  the  Sun  and  Moon.  The 
periods  are  close  to  diurnal  because  the  Sun  and  Moon  rise  and  set  once  a  day. 

For  many  applications  nutations  are  a  nuisance  since  if  they  are  not  adequately  removed  they  can 
degrade  solutions  for  other  parameters.  However,  recent  observational  results  from  Very-Long- 
Baseline-Interferometry  (VLBI)  have  demonstrated  that  nutations  can  be  important  in  their  own  right  in 
providing  a  probe  of  the  Earth’s  interior.  In  this  chapter,  we  review  nutation  theories  and,  in  an  appen¬ 
dix,  we  extend  those  theories  to  include  an  Earth  with  a  non-hydrostatic  equilibrium  state.  We  discuss 
the  geophysical  implications  of  the  observational  results. 


A  REVIEW  OF  THE  THEORY 


The  goal  of  a  nutation  theory  is  to  estimate  the  nutation  amplitude  as  a  function  of  frequency. 
The  amplitude  depends  on  frequency  because  the  forcing  from  the  Sun  and  Moon  depends  on  fre¬ 
quency,  and  because  the  Earth’s  response  to  forcing  is  different  at  different  frequencies.  Finding  the 
forcing  as  a  function  of  frequency  (in  fact,  finding  the  frequencies  themselves)  is  a  celestial  mechanics 
problem,  and  requires  accurate  solutions  for  the  orbital  motion  of  the  Moon  and  Earth  [Kinoshita, 
1977J.  Modelling  the  Earth’s  response  as  a  function  of  frequency  is  a  geophysical  problem,  and  it  is 
that  problem  that  will  be  addressed  in  this  chapter. 

The  gravitational  potential  energy  from  the  Sun  and  Moon  can  be  expanded  in  Earth-based  coor¬ 
dinates  as  a  sum  of  complex  spherical  harmonics  (T").  Each  coefficient  in  that  sum  can  be  further 
expanded  as  a  sum  of  terms  which  vary  harmonically  with  time.  Only  Y 2  terms  contribute  to  the  nuta¬ 
tional  motion.  We  write  an  individual  Yj  term  as 

V  =/((o)e"“r2K21(8,X)  (1) 


where  to  is  the  frequency,  0  and  k  are  the  co-latitude  and  eastward  longitude,  /  (to)  is  a  scalar  ampli¬ 
tude,  and  Y 2  is  normalized  so  that  the  integral  of  IY2  I2  over  the  unit  sphere  is  1.  The  values  of  to  and 
/  (to)  are  determined  by  the  orbital  motion. 

Under  very  general  conditions,  the  response  of  the  Earth  to  the  potential  (1)  can  be  expanded  as  a 
sum  of  the  Earth’s  normal  mode  eigenfunctions  [Wahr,  1981a],  so  that  the  nutation  amplitude,  £  is 


C  =  e,<0,Z 


B: 


-w-co, 


(2) 


Here,  the  sum  over  i  is  over  all  normal  modes,  to,  is  the  eigenfrequency  of  the  i’th  mode,  and  5, 
depends  on  the  i’th  eigenfunction  and  (linearly)  on  / .  One  implication  of  (2)  is  that  the  amplitude  is 
large  if  the  forcing  frequency,  to,  is  close  to  an  eigenfrequency. 

Although  the  sum  in  (2)  is,  in  principle,  over  every  one  of  the  Earth’s  normal  modes,  almost  all 
of  the  important  contributions  come  from  just  two  modes,  both  with  frequencies  in  the  diurnal  band. 
One  of  these  modes,  called  the  tilt-over-modc  (TOM),  is  simply  a  tipping  of  the  Earth  in  space  with  no 
associated  deformation.  The  motion  is  exactly  equivalent  to  lipping  the  coordinate  system  in  the  oppo¬ 
site  direction.  The  period  of  the  TOM  is  infinite  as  seen  from  non-rotating  inertial  space,  and  so  is 
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exactly  one  sidereal  day  as  seen  from  a  sidereally  rotating  system. 

The  more  interesting  of  the  two  diurnal  modes  contributing  to  (2)  is  the  free  core  nutation  (FCN). 
This  mode  involves  tipping  the  mantle  and  fluid  core  in  opposite  directions.  Because  the  core-mantle 
boundary  is  not  exactly  spherical,  the  mantle  and  core  push  against  each  other  when  they  tip,  and  the 
resulting  pressure  acts  to  restore  the  core  and  mantle  to  their  untipped  state.  There  are  also  gravita¬ 
tional  restoring  forces  due  to  the  interaction  between  the  aspherical  density  distributions  of  the  core  and 
mantle.  The  result  is  a  free  periodic  motion,  the  FCN,  with  a  frequency  equal  to  one  cycle  per  day  plus 
a  small  term  dependent  on  the  strength  of  the  core-mantle  coupling.  That  coupling  depends  on  the 
aspherical  shape  and  density  distribution  of  the  core,  as  described  below.  For  an  Earth  where  the  shape 
and  density  are  assumed  to  be  consistent  with  a  state  of  hydrostatic  pre-stress,  the  frequency  is 
torav  =  1+1/460  c/d  [Wahr,  1981b]. 

The  TOM  contributions  to  the  nutation  amplitudes  in  (2)  are  much  larger  than  the  FCN  contribu¬ 
tions.  The  reason  is  that  for  the  TOM,  the  factor  B,  in  (2)  depends  on  the  total  torque  on  the  Earth 
from  the  Sun  and  Moon  at  the  frequency  to.  For  the  FCN,  B,  is  sensitive,  instead,  to  the  difference 
between  the  torques  per  unit  moment  of  inertia,  acting  on  the  core  and  mantle  (and  to  deformation 
terms  of  about  the  same  order).  The  Sun  and  Moon  do  not  provide  nearly  as  large  a  differential  torque 
as  they  do  a  total  torque. 

The  TOM  contributions  are  not  particularly  interesting.  The  TOM  resonance  merely  reflects  the 
fact  that  the  longer  (  as  seen  from  inertial  space)  you  torque  an  object,  the  more  it  tips.  The  FCN  con¬ 
tributions,  though,  are  of  geophysical  importance  since  they  are  sensitive  to  the  poorly  known  shape  and 
asplierical  density  distribution  of  the  core.  To  consider  those  contributions  in  detail,  it  is  usual  to 
remove  the  TOM  contribution  from  the  sum  in  (2),  and  then  to  divide  the  remainder  by  the  TOM  con¬ 
tribution.  This  last  division  removes  the  scalar  amplitude,  /  (to),  from  the  remainder,  so  that  the  results 
reflect  only  the  geophysically  interesting  part  of  the  signal.  This  normalized  "admittance"  is  shown  in 
Figures  III.l  and  III.2  as  a  function  of  frequency.  The  theoretical  results  (the  solid  line)  use  eigenfre- 
quencies  and  eigenfunctions  computed  for  a  hydrostatically  pre-stressed  Earth.  Although  the  results 
shown  are  computed  using  several  modes  in  the  sum  (2)  [Wahr,  1981b],  almost  all  of  the  contributions 
come  from  the  FCN.  The  resonance  at  the  FCN  cigenfrcquency  of  1+1/460  c/d  is  clearly  evident. 

Also  shown  in  Figures  III.l  and  III.2  are  recent  VLBI  observational  results  [Herring  et  al.,  1986] 
for  the  nutation  amplitudes  at  a  few  important  frequencies.  The  TOM  resonance  has  been  removed 
from  the  observational  results  and  then  divided  into  the  remainder,  in  order  to  compare  with  the  theoret¬ 
ical  results  (the  observations  have  also  been  corrected  for  the  effects  of  oceans  as  described  below). 
The  vertical  bars  on  the  observational  results  reflect  the  published  errors. 

The  agreement  between  the  theory  and  the  observations  is,  in  general,  good.  The  disagreement  at 
the  prograde  fortnightly  frequency  is  observational  ly  significant  but,  in  absolute  terms,  is  only  a  few 
tenths  of  a  milli-arcsecond.  The  agreement  at  the  prograde  semi-annual  frequency  looks  from  Figure 

111.1  to  be  reasonably  good,  but,  as  we  shall  see  below,  there  could  be  a  disagreement  of  up  to  1  milli- 
arcsecond  after  correcting  for  mantle  anelasticity  and  non-hydrostalic  core  structure. 

First,  however,  there  is  an  even  larger  discrepancy  at  the  retrograde  annual  frequency  (see  the 
enlarged  comparison  in  Figure  111.2).  The  observed  annual  admittance  lies  well  above  the  theoretical 
result,  and  the  difference  in  absolute  terms  is  about  2  milli-arcseconds.  The  results  in  Figures  III.l  and 

111.2  suggest  that  the  FCN  eigenfrequency  should  be  larger  than  the  value  of  1+1/460  c/d  predicted 
using  the  hydrostatic  assumption.  In  fact,  Gwinn,  ct  al.  [1986]  used  the  observational  results  to  con¬ 
clude  that  (arcs  =  1+1/433  c/d.  (Although  the  annual  discrepancy  could  also  be  resolved  by  adjusting 
the  FCN  value  of  fi,  in  (2),  that  adjustment  would  cause  substantial  discrepancies  al  other  frequencies.) 

What  does  this  increase  in  eigenfrequency  imply  about  the  Earth?  The  theoretical  results  shown 
in  Figures  III.l  and  III.2  do  not  include  the  effects  of  oceans  (although  the  observational  results  have 
been  corrected  for  the  oceans),  of  mantle  anelasticity,  or  of  non-hydrostalic  pre-stress  and  structure. 
These  effects  are  discussed  below. 
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The  Oceans 


Oceans  allect  nutation  amplitudes  through  surface  loading.  The  Sun  and  Moon  cause  diurnal 
tides  in  the  oceans  at  exactly  the  nutation  periods.  Those  ocean  tides  load  the  Earth  and  cause  further 
nutational  motion.  Thus,  the  effects  of  oceans  can  be  perceived  as  modifying  the  driving  force  for 
nutations.  In  that  case,  the  force  can  no  longer  be  written  in  terms  of  a  Y 2  potential,  as  in  (1).  Oce¬ 
anic  corrections  require  some  understanding  of  the  loading  force,  and  that  requires,  at  the  very  least, 
reliable  ocean  tide  models.  Wahr  and  Sasao  [1981]  used  diurnal  tide  models  to  estimate  the  oceanic 
corrections,  and  their  results  have  been  removed  from  the  VLBI  observations  to  give  the  results  shown 
in  Figures  III.  I  and  III.2. 

Mantle  Anelasticity 


Mantle  anelasticity  and  non-hydrostatic  structure  affect  the  nutations  by  modifying  the  Earth’s 
response  to  external  forcing,  rather  than  by  contributing  to  the  forcing  itself.  To  understand  their  con¬ 
tributions,  it  is  necessary  to  describe  nutation  models  in  more  detail. 

Existing  nutation  models  fall  into  two  categories,  referred  to  here  as  numerical  and  semi- 
analytical.  Both  types  of  models  involve  the  solution  of  the  same  infinite  set  of  coupled  ordinary 
differential  equations.  And,  both  types  of  models  derive  approximate  solutions  by  truncating  the  equa¬ 
tions.  In  the  numerical  method  [Wahr,  1981b;  see  also  Smith,  1977]  the  truncation  is  less  severe  and 
more  terms  are  kept  than  in  the  semi-analytical  method  (see,  for  example,  Jeffreys  and  Vicente  [1957]; 
Molodensky  [1961];  Sasao  et  al.  [1980]).  Thus,  the  advantage  of  the  numerical  method  is  that  it  is  apt 
to  be  more  accurate.  The  advantages  of  the  semi-analytical  methods  are  (1)  they  are  easier  to  imple¬ 
ment;  and  (2)  the  results  are  more  readily  understood  in  terms  of  the  Earth’s  physical  parameters. 
Furthermore,  the  results  from  the  semi-analytical  models  appear  to  agree  well  with  those  from  the 
numerical  method. 

The  semi-analytical  results  of  Sasao  et  al.  [1 980}  for  BFCn  and  coFCN  are: 


Bfcn  =  -VT5/2S 


<0FCN  - 


Q 


(3) 

(4) 


where  A,  Af,  Am  are  the  equatorial  moments  of  inertia  of  the  Earth,  core,  and  mantle;  e  =  {C-A  )IA 
and  eF  =  (CF  -Af  )IAF  are  the  dynamical  elliplicities  of  the  Earth  and  core  (C  and  Cy  are  the  polar 
moments  of  inertia  of  the  Earth  and  core);  and  y  and  P  represent  the  effects  of  deformation  within  the 
mantle  and  core  (y<e  and  P<«/).  The  results  shown  in  (3)  and  (4)  were  derived  by  Sasao,  et  al.  [1980] 
assuming  the  inner  core  is  fluid,  and  that  the  Earth's  equilibrium  state  is  one  of  hydrostatic  pre-stress. 

Mantle  anelasticity  affects  these  results  by  modifying  y  and  p.  Those  parameters  become  com¬ 
plex,  leading  to  phase  lags  in  the  nutations.  Even  larger  are  the  effects  on  the  real  part  of  those  param¬ 
eters.  The  main  source  of  uncertainty  when  modelling  the  anelastic  corrections  is  the  uncertainty  in 
mantle  Q  at  diurnal  periods.  Seismic  information  is  pertinent  to  much  shorter  periods.  On  the  other 
hand,  this  suggests  that  perhaps  the  nutation  results  could  be  used  to  leant  about  mantle  Q  in  this  fre¬ 
quency  regime. 

Wahr  and  Bergen  [1986]  and  Dehant  [1988]  used  a  variety  of  anelastic  models  and  assumptions 
to  estimate  the  effects  on  nutation  amplitudes.  They  found  that,  in  general,  the  contributions  to  BFCN  in 
(2)  were  more  important  than  the  contributions  to  c &Fcn-  The  model  results  shown  in  Figures  III.  1  and 
III.2  are  corrected  for  Wahr  and  Bergen’s  anelastic  estimates  and  the  results  are  shown  in  Figures  III. 3 
and  III.4.  Also  shown,  in  Figure  III. 5,  is  a  comparison  between  the  oul-of-phase  components  inferred 
from  the  VLBI  results  (corrected  for  the  effects  of  oceans)  and  the  predicted  results  due  to  anelasticity. 
For  a  dissipationless  Earth,  the  out-of-phase  components  would  be  zero.  The  vertical  bars  associated 
with  the  anelastic  results  in  these  figures  reflect  the  uncertainty  in  the  value  of  mantle  Q  at  diurnal 
periods.  The  figures  suggest  that  the  present  VLBI  results  are,  in  principle,  nearly  accurate  enough  to 
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discriminate  between  various  mantle  Q  results.  However,  the  effects  of  anelasticity  do  not  resolve  the 
disagreements  between  the  observations  and  the  theory.  In  fact,  they  tend  to  make  the  agreement  for 
the  in-phase  amplitudes  (Figs.  3  and  4)  worse.  Thus,  before  the  nutation  results  can  be  used  to  leam 
about  anelasticity,  it  will  first  be  necessary  to  understand  the  reasons  for  the  large  discrepancy. 


Non-hydrostatic  Structure 


The  most  likely  explanation  for  the  difference  between  the  observed  and  theoretical  cofCW  is  that 
the  Earth  is  not  hydrostatically  pre-stressed.  In  the  Appendix,  we  extend  the  semi-analytical  model  of 
Sasao  et  al.  to  include  an  Earth  which  is  not  hydrostatically  pre-stressed  (although  we,  too,  assume  the 
inner  core  is  fluid).  In  particular,  the  core-mantle  boundary  and  density  distribution  within  the  Earth 
can  have  an  arbitrary  spherical  harmonic  expansion.  We  find,  though,  that  BFCN  and  u>FCN  are  still 
given  by  equations  (3)  and  (4).  Thus,  the  effect  of  the  non-hydrostatic  pre-stress  on  (0FCN  can  be  com¬ 
puted  by  estimating  the  non -hydrostatic  contributions  to  ef  in  (4).  Since  e;  depends  only  on  the  T® 
spherical  harmonic  components  of  the  core-mantle  boundary  shape  and  of  the  core’s  internal  density, 
the  observational  results  for  cofc/v  constrain  only  those  particular  spherical  harmonic  coefficients. 


The  VLBI  results  for  (ssFCn  suggest  that  et  is  about  5%  larger  than  the  hydrostatic  value  [Gwinn 
et  al.,  1986].  Suppose  the  radius  of  the  core-mantle  boundary  as  a  function  of  8  and  X  is 


r(0,X)  =  a 


1  +  L5>r>mx) 

1=1  *=-/ 


(5) 


where  a  is  the  mean  core  radius,  and  the  Y,m  are  normalized  so  that  the  integral  of  I >7” 1 2  over  the  unit 
sphere  is  1.  We  have  found  numerically  [J.  Wahr  and  D.  de  Vries,  unpublished  manuscript],  that  unless 
there  is  a  thin,  low  density  boundary  layer  at  the  top  of  the  core,  ef  is  approximately  equal  to  the  ellip- 
ticity  of  the  core  mantle  boundary,  ec,  where  a®  =  -(2/3 )ec.  (If  there  is  a  low  density  boundary  layer, 
then  ef  could  also  depend  critically  on  the  shapes  of  the  equi-density  surfaces  inside  the  core.)  Thus,  in 
the  absence  of  such  a  boundary  layer,  the  nutation  results  constrain  the  K®  component  of  the  boundary. 
They  imply,  in  fact,  that  tc  is  about  5%  larger  than  the  hydrostatic  value,  and  suggest  a  non-hydrostatic 
topography  on  the  boundary  of  about  1/2  km. 

When  the  value  usFCN  =  1  +■  1/433  c/d  is  used  in  (2)  in  place  of  the  hydrostatic  result  1  +  1/460 
c/d,  large  differences  occur  at  other  frequencies,  as  well.  For  example  [Gwinn,  et  al.,  1986],  the 
theoretical  retrograde  18.6-year  amplitude  is  decreased  relative  to  the  IAU  adopted  value  by  about  2 
milli-arcsecond.  This  amplitude  is  also  decreased  by  about  1  milli-arcsecond  due  to  the  oceans  and 
about  .5  milli-arcsecond  by  anelasticity. 


The  theoretical  prograde  six-month  term  is  increased  by  about  .4  milli-arcseconds  due  to  the 
change  in  u>FCN.  In  fact,  the  VLBI  result  for  this  term  disagrees  with  the  theoretical  amplitude  by 
about  1  milli-arcsecond,  after  correcting  for  the  effects  of  the  oceans  and  anelasticity,  and  using  aFCN  = 
1  +  1/433  c/d.  Similarly,  there  is  a  discrepancy  for  the  prograde  fortnightly  term  of  about  .3  milli- 
arcseconds  between  the  VLBI  results  and  the  adjusted  theory.  It  is  not  clear,  at  present,  what  is  respon¬ 
sible  for  the  fortnightly  and  semi-annual  discrepancies. 


CHAPTER  III  REFERENCES 
Dahlen,  F.A.,  Geophys.  J.  Roy.  Astr.  Soc..  32,  203,  1973. 

Dehant,  V.,  in  Earth  Rotation  and  Reference  Frames,  eds.  Babcock,  A.,  and  G.A.  Wilkins,  1988. 
Gwinn,  C.R.,  T.A.  Herring,  and  1. 1.  Shapiro,  J.  Geophys.  Res.,  91,  4755,  1986. 

Herring,  T.A.,  C.R.  Gwinn,  and  I.l.  Shapiro,  J.  Geophys.  Res.,  91,  4745,  1986. 


22 


Jeffreys,  H„  and  R.O.  Vicente,  Mon.  Not.  R.  Astr.  Soc.  117,  162,  1957. 

Kinoshita,  H„  Cel.  Mech.,  15,  277,  1977. 

Lieske,  J.,  T.  Lederle,  W.  Fricke,  and  B.  Morando,  Astron.  Astrophys.,  58,  1,  1977. 

Molodensky,  M.S.,  Commun.  Obs.  R.  Belg.,  288 , 25,  1961. 

Saito,  M„  J.  Phys.  Earth.  22,  123,  1974. 

Sasao,  T.,  S.  Okubo,  and  M.  Saito,  in  Proceedings  of  IAU  Symposium  no.  78  "Nutation  and  the  Earth's 
Rotation",  eds.  Fedorov,  E.,  M.  Smith,  and  P.  Bender,  D.  Reidel  Publishing  Co.,  Dordrecht,  1980. 

Smith,  M.L.,  Geophys.  J.  Roy.  Astr.  Soc.,  50,  103,  1977. 

Smith,  M.L.,  and  F.A.  Dahlen,  Geophys.  J.  Roy.  Astr.  Soc.,  64,  223,  1981. 

Wahr,  J.M.,  Geophys.  J.  Roy.  Astr.  Soc.,  64,  651,  1981a. 

Wahr,  J.M.,  Geophys.  J.  Roy.  Astr.  Soc.,  64,  705,  1981b. 

Wahr,  J.M.  and  T.  Sasao,  Geophys.  J.  Roy.  Astr.  Soc.,  64,  747,  1981. 

Wahr.  J.,  and  Z.  Bergen,  Geophys.  J.  Roy.  Astr.  Soc.,  87,  633,  1986. 

Wahr,  J.,  and  D.  dc  Vries,  Geophys.  J.  Roy.  Astr.  Soc.,  1989  (submitted). 

Woodhouse,  J.H.,  and  F.A.  Dahlen,  Geophys.  J.  Roy.  Astr.  Soc.,  53,  335,  1978. 


CHAPTER  III  APPENDIX:  Nutation  Amplitudes  for  a 
Non-hydrosiaiicaliy  Pre-stressed  Earth 

Previous  nutation  models  have  assumed  the  Earth  is  hydrostatically  pre-stressed.  Among  the  implica¬ 
tions  of  this  assumption  is  that  the  core-mantle  boundary  and  all  surfaces  of  constant  density  within  the 
core  and  mantle  arc  elliptical.  Thus,  the  departure  of  those  surfaces  from  spherical  symmetry  is  propor¬ 
tional  to  the  single  spherical  harmonic  Y°(Q,\). 

It  is  quite  likely,  though,  that  surfaces  of  constant  density  in  the  Earth  and  surfaces  >f  discontinuity, 
such  as  the  core-mantle  boundary,  have  non -hydrostatic  shapes.  The  Y°  components  of  these  surfaces 
are  apt  to  differ  from  the  hydrostatic  values,  and  there  are  apt  to  be  components  with  other  Y”  angular 
dependence.  Since  the  FCN  is  critically  dependent  on  pressure  coupling  acting  across  the  core-mantle 
boundary  (and,  to  a  lesser  extent,  gravitational  coupling  between  the  core  and  mantle),  the  FCN  contri¬ 
bution  to  the  eigenfunction  expansion  (2)  could  be  affected. 

In  this  appendix,  the  FCN  contribution  is  modelled  without  assuming  hydrostatic  equilibrium  in  the 
mantle.  Instead,  the  mantle  is  allowed  to  possess  a  small,  arbitrary  deviatoric  pre-stress,  and  so  surfaces 
of  constant  material  properties  and  surfaces  of  discontinuity  (such  as  the  core/mantle  boundary)  within 
the  mantle  can  have  arbitrary  shape.  The  core,  however,  is  not  permitted  to  have  a  deviatoric  pre-stress 
(the  core,  after  all,  is  a  fluid  even  at  seismic  periods).  As  a  result,  the  density  distribution  in  the  core  is 
uniquely  determined  by  the  shape  of  the  core-mantle  boundary  and  by  the  gravitational  potential  from 
the  mantle  acting  on  the  core.  Thus,  the  constant  density  surfaces  in  the  core  need  not  coincide  with 
the  surfaces  predicted  for  an  everywhere  hydrostatically  pre-stressed  Earth.  The  derivation,  below, 
extends  the  semi-analytical  model  of  Sasao,  el  al.  [1980]  to  include  these  non-hydrostatic  modifications, 
and  to  more  fully  justify  a  number  of  the  necessary  approximations.  As  in  the  model  of  Sasao,  et  al. 
[1980],  we  assume  the  inner  core  is  fluid. 

Consider  a  reference  frame  fixed  to  the  mantle  and  rotating  with  it.  For  nutational  motion  at  fre¬ 
quency  to,  the  rotation  vector  of  the  frame  is  defined  as  +  Qm0(x+iJ)e‘a>',  where  the  nutational 
motion  of  the  mantle  is  described  by  m0,  and  where  m0«l.  For  the  forced  nutations,  the  ratio 
(ci)-fJ)/Q  is  small.  We  will  consider  this  ratio  to  be  first  order  in  E,  where  e  is  some  measure  of  the 
departure  of  the  Earth  from  spherical  symmetry.  This  assumption  is  reasonable  for  the  annual  nutations 
(the  terms  with  the  largest  discrepancy  between  theory  and  observation),  where  l(o)-nyQ  1=1/365.  It  is 
less  accurate  for  the  fortnightly  nutations,  where  l(o>-Q)/ftl=l/13.7. 

Denote  the  luni-solar  tidal  force  by  -pW,  where  V  is  the  luni-solar  tidal  potential  given  by  equation 
(1)  in  the  text.  The  equations  of  motion  in  the  hydrostatically  pre-stressed  core  arc  [Sasao,  cl  al.. 
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1980]: 


pI-co2s+2'0)Qxs]  =  -pV(0,H4'}-VP£-piV<t>-ft2pmo[(x+jy)z(£}-a)) 

+z(x  +iy  )(£2-hd)]  (A  1 ) 

V24>,  =  4kG  p, 

where  s  is  the  displacement  field  in  the  core;  0lt  PE,  and  P)  arc  the  Eulerian  perturbations  in  gravita¬ 
tional  potential  energy,  pressure,  and  density,  respectively;  is  the  mean  rotation  vector  of  the 

Earth;  and  p  and  <t>  are  the  initial  density  and  total  potential  energy  (gravitational  plus  centrifugal)  and 
can  have  arbitrary  angular  dependence,  except  that  surfaces  of  constant  p  must  be  surfaces  of  constant 
d>  in  the  core.  There  are,  also,  equations  relating  P)  and  PE  to  s,  which  are  not  shown  here.  And, 
there  are  similar,  although  more  complicated,  differential  equations  describing  displacements  in  the 
mantle,  which  arc  also  omitted  here.  The  mantle  equations  must  explicitly  include  any  assumed  non¬ 
hydrostatic  pre-stress  (see  Woodhousc  and  Dahlen  (1978]  for  the  complete  mantle  and  core  equations 
•md  for  the  boundary  conditions). 

For  nutational  motion,  the  displacement  field  in  the  core  has  the  form 

s(r)  =  f20o(x+i  y)xr  +  s^r)  (A2) 

where  60  represents  the  mean  nutational  motion  relative  to  the  mantle,  and  sd,  which  represents  the 
deformation  of  the  core,  is  first  order  in  e  compared  with  the  90  term  (see,  for  example,  Sasao,  et  al. 
[1980]).  T-.-ing  (x-iy)  [rx(equation(Al))]  and  integrating  through  the  core,  and  using  (A2)  for  s,  gives 
an  angular  momentum  equation  for  the  core: 

2iQ2Afm0  +  20iQ2AfQo~^- 

=  (: x-iy )•  |  rx  -pV^-t-l/^V/^-p,  Vd>  .  (A3) 

core  J 

Here,  Af  is  the  core’s  principal  moment  of  inertia  in  the  equatorial  plane;  and  cL=c{i~ich ,  where  the 
c(t  are  perturbations  in  the  core’s  inertia  tensor  due  both  to  the  deformational  displacements,  sd,  and  to 
the  rotational  motion,  0O  (the  components  of  the  inertia  tensor  change  when  the  aspherical  core  is 
rotated  with  respect  to  the  coordinate  system).  The  terms  on  the  right  hand  side  of  (A3)  include  the 
pressure  and  gravitational  torques  on  the  core  from  the  mantle,  and  the  gravitational  torque  on  the  core 
from  the  luni-solar  tidal  force. 

All  terms  in  (A3)  are  second  order  in  e  relative  to  90.  As  an  example,  <J>)  has  contributions  from  the 
deformation,  sd ,  and  from  the  rotational  motion,  described  by  80.  The  contributions  to  <J>i  from  the 
deformation  are  first  order,  because  sd  is  first  order.  The  contributions  from  0O  are  caused  by  rotating 
the  Earth’s  unperturbed  gravity  field  through  the  angle  90.  These  rotational  contributions  are  also  first 
order,  because  there  is  no  change  in  gravity  if  the  unperturbed  field  is  spherically  symmetri  .  Thus,  <]>i 
is  a  first  order  quantity.  Furthermore,  the  integral  of  prxV0j  is  0  no  matter  what  4>i  is,  if  both  the  core 
shape  and  the  internal  density  in  the  core  are  spherically  symmetric.  Thus,  J  rxpVtjq  is  second  order. 

core 

As  a  second  example,  fioefl2  is  a  first  order  quantity  (the  unperturbed  centrifugal  force  is  order  e 
times  the  unperturbed  gravitational  force),  as  is  (f2-cu)/£2.  Thus,  the  0O  term  in  (A3)  is  second  order. 
And,  as  a  third  example,  m0  is  first  order  compared  with  Q0O  (m0  represents  motion  of  the  mantle  rota¬ 
tion  axis  as  seen  from  the  mantle  and  is  of  order  (f2-e))0o)  so  that  Q2me  is  second  order.  Finally,  for 
completeness,  we  note  here  that  V,  PE,  pj,  <t>,  and  cL  arc  all  first  order  rclativp  to  0O. 

A  comparable  angular  momentum  equation  for  the  entire  Earth,  after  dropping  all  terms  third  order  or 
smaller,  is 

2 i  Q2Am0  =  2(<o-f2)A/  i  (C-A)  f  (to)  (A4) 

where  C  and  A  are  the  principal  moments  of  inertia  for  the  entire  Earth,  /  (to)  is  the  scalar  amplitude 
of  the  tidal  potential  (see  equation  (1)  in  the  text),  and  the  last  term  on  the  right  hand  side  of  (A4) 
represents  the  luni-solar  torque  on  the  Earth. 
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Next,  cL  and  the  right  hand  side  of  (A3)  can  be  related  to  0O.  Use  (A2)  in  (Al)  and  separate  the 
resulting  vector  equations  into  spheroidal  and  toroidal  scalar  equations.  Similarly,  separate  the  mantle 
differential  equations  and  the  boundary  conditions  into  toroidal  and  spheroidal  scalar  equations.  Then, 
consider  only  the  spheroidal  equations,  ignore  all  terms  in  these  equations  that  are  second  order  or 
smaller  in  e,  and  solve  the  entire  system  on  a  computer.  In  the  core,  for  example,  the  first  order 
spheroidal  equations  derived  from  (A  I)  arc  the  scalar  components  of 

p0V«>i+V>+V/’£+p1V<J>()=  [/0on3V8S7l5]  p0V(r2y2')  (A5) 

where  p0  and  <t>0  are  the  spherically  symmetric  parts  of  p  and  <t>.  (One  of  the  consequences  of  the  trun¬ 
cation  to  first  order  is  that  m0  does  not  appear  in  (A5).) 

The  first  order  differentia!  equations  in  the  mantle  corresponding  to  (A5)  in  the  core,  are  the  usual  set 
of  spheroidal  equations  describing  tidal  deformation  of  a  spherical,  non-rotating,  static,  and  hydrostati¬ 
cally  pre -stressed  mantle  (see,  for  example,  Saito  [1974]).  Although,  in  principle,  there  is  an  apparent 
spheroidal  force  in  the  mantle  which  depends  on  m0,  that  force  is  second  order  and  so  can  be  ignored 
in  the  first  order  deformation  equations.  All  first  order  boundary  conditions  within  the  core  and  mantle 
and  at  the  outer  surface  are  also  equivalent  to  the  boundary  conditions  for  a  spherical,  non-rotating, 
hydrostatically  pre-stressed  Earth. 

In  effect,  then,  all  deformation  terms  can  be  computed  by  solving  the  static  equations  of  motion  for  a 
spherical,  non-rotating,  hydrostatically  pre-stressed  Earth,  subject  to  an  apparent  force  proportional  to 
p0V(r2y2 ),  but  with  different  proportionality  constants  in  the  core  and  mantle.  In  the  core,  the 
apparent  force  is  proportional  to  0O  and  /  (to)  (see  (1)  for  a  definition  of  /  (to)).  In  the  mantle,  the 
apparent  force  is  proportional  only  to  /  (to).  Neither  the  non-hydrostatic  pre-stress  in  the  mantle  nor 
any  of  the  Earth's  asphcrical  structure  enters  explicitly  into  any  of  the  first  order  deformation  equations. 
Their  effects  are  included  only  through  the  integrals  on  the  right  hand  side  of  (A3). 

Because  the  apparent  force  in  the  core  is  proportional  to  p0V(r2y2 )  and  because  the  first  order  defor¬ 
mation  equations  are  spherically  symmetric,  PE ,  and  <j>|  have  K]  (0.X)  angular  dependence.  Using 
this  angular  dependence,  (A5)  yields  directly 

P£(r,0,X)  =  p0(r ) (t Q30on'8tc7T5-/ (to)] r 2Y±  (0,X)-<l>i(r  ,0,X>]  (A6) 

Pi(r  ,0,X)  =  [<»,(/■  ,0.X)+(/  (0))-/ f230o«5]  r2h  >  (0A)]  . 


Using  (A6)  in  (A3)  and  doing  the  integrals  gives  a  result  accurate  to  second  order  of 

,  ,  iciSl2  , 

2iQ2m0+2n\(n-(o}+ - =  -2Q\e,  (A7) 

Af 

where  ef=(Cf-Af)IAf  is  the  dynamical  eliipticily  of  the  core  (Cj-  is  the  polar  moment  of  inertia  of  the 
core). 

Define  the  dimensionless,  real  parameters  (3  and  y  so  that 

cL  =  (32 A;  i  QQo-yAf  sfl572n  /  (to )/&  .  (A8) 


(3  and  y  can  be  determined  by  solving  the  deformation  equations  on  a  computer.  Using  (A4)  and  (A8) 
in  (A7)  and  solving  for  the  core  rotation  angle  0O,  gives 


where  e=(C-A)/A  is,  to  lowest  order,  the  dynamical  ellipticity  of  the  Earth,  and  Am  is  the  principal 
moment  of  inertia  of  the  mantle. 


The  nutation  amplitude  observed  at  the  Earth’s  outer  surface  is  (see  Sasao  and  Wahr  (1981,  eq.  3.20]) 
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(A10) 


C  = 


£1 

-x — m0  . 
ft— <o 


Using  (A4)  to  relate  mQ  to  90,  and  (A9)  to  relate  0O  to  /  (to),  (A10)  reduces  to 


c= 


e 

£2-0) 

' 

CD- 

I+JL^-P) 

Vl5 . 


(All) 


The  e£2/(£2-co)  term  in  (All)  represents  the  TOM  resonance.  The  other  term  in  (All)  is  the  FCN 
resonance,  and  can  be  written  as  BFCNl(t»-u>FCfj),  where  BFCN  and  (aFCN  are  given  by  equations  (3)  and 
(4)  in  the  text. 

Although  this  result  was  derived  here  without  assuming  a  hydrostatically  pre-stressed  mantle,  it  is 
identical  in  form  to  the  hydrostatic  result  The  dynamical  ellipticity,  ,  depends  on  the  T®  component 
of  the  core-mantle  boundary  shape,  and  on  the  T®  terms  in  the  density  stratification  inside  the  core. 
The  dependence  of  (£Fcn  on  the  internal  density  stratification  is  due  to  the  effects  of  gravitational 
torques  between  the  core  and  mantle,  represented  by  the  pVtpj  and  pi  V$  terms  on  the  right  hand  side  of 
(A3).  The  dependence  on  the  boundary  structure  is  due  to  pressure  torques  at  the  core/mantle  boun¬ 
dary,  represented  by  the  VPE  term  in  (A5).  There  is  no  dependence,  to  this  order  of  approximation,  on 
any  other  Y"  terms  in  the  aspherical  structure.  Similarly,  e  depends  on  the  K®  density  structure 
throughout  the  Earth,  and  is  well  determined  from  independent  observations  of  the  Earth’s  precession 
(see,  for  example,  Lieske  et  al.  [1977]).  The  factors  P  and  y  represent  the  effects  of  deformation  and 
are  insensitive,  to  this  order,  to  aspherical  structure.  Sasao  et  al.  [1980]  found  that  p  is  about  25%  of 

the  hydrostatic  value  of  ef.  For  a  hydrostatically  pre-stressed  Earth,  coravs(l+^^)  cycles  per  day. 
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CHAPTER  III  FIGURE  CAPTIONS 


Figure  III.  1.  A  comparison  between  theoretical  nutation  results  from  Wahr  (1981b)  (solid  line)  and  VLBI 
observational  results  from  Herring,  et  al.  (1986)  (vertical  bars).  The  results  are  normalized,  as  described  in 
the  text,  to  form  admittances.  The  lengths  of  the  vertical  bars  represent  the  observational  errors. 

Figure  III.2.  An  expanded  view  of  the  1+1/365  term  in  Figure  III.l.  There  is  significant  disagreement 
between  theory  and  observation.  The  discrepancy  suggests  the  FCN  eigenfrequency  should  be  larger  than 
the  theoretical  value. 

Figure  III.3.  The  theoretical  nutation  admittances  for  an  anelastic  Earth  from  Wahr  and  Bergen  (1986). 
The  vertical  bars  for  the  anelastic  results  reflect  the  uncertainty  in  mantle  Q  at  diurnal  periods.  Also  shown 
are  the  theoretical  results  for  an  elastic  Earth  (solid  line)  and  the  VLBI  results. 

Figure  III.4.  The  anelastic  admittance  and  observational  result  for  the  I+I/365  nutation.  Anelasticity  wor¬ 
sens  the  agreement  between  observation  and  theory. 

Figure  III.5.  The  out-of-phase  components  for  the  nutations.  Both  the  VLBI  results  (solid  lines)  corrected 
for  the  oceans  and  the  anelastic  contributions  (dashed  lines)  are  shown.  The  vertical  bars  for  the  anelastic 
results  reflect  uncertainty  in  mande  Q  at  diurnal  periods.  The  out-of-phase  components  should  be  zero  for 
a  dissipationless  Earth. 
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Figure  M.2 
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CHAPTER  IV 


THE  POSSIBILITY  OF  LATERAL  STRUCTURE  INSIDE  THE  CORE  AND  ITS 
IMPLICATIONS  FOR  NUTATION  AND  EARTH  TIDE  OBSERVATIONS 


SUMMARY 

Wc  show  that  there  can  be  lateral  structure  inside  the  fluid  core  caused  by  gravitational  forcing  from 
the  mantle,  from  the  inner  core,  or  from  topography  on  the  corc/mantle  or  inner  core/outer  core  boun¬ 
daries.  We  describe  a  method  for  calculating  the  internal  structure,  given  knowledge  of  the  forcing. 
We  discuss  the  possible  effects  of  this  structure  on  results  for  the  core/mantle  boundary  topography 
inferred  from  observations  of  the  earth’s  forced  nutations  and  diurnal  earth  tides.  We  consider  the  pos¬ 
sible  implications  of  a  thin,  low  density  fluid  layer  at  the  top  of  the  core. 

1  INTRODUCTION 

There  has  been  considerable,  recent  effort  in  the  geophysics  community  to  understand  more  about  the 
structure  and  the  dynamical  and  material  properties  at  the  core/mantle  boundary.  Constraints  on  the 
boundary  shape  and  on  the  near-boundary  structure  come  from  numerous  sources,  including  seismic 
imaging  with  body  waves  [eg.  Gudmundsson,  et  al„  1986;  Creager  and  Jordan,  1986;  Morelli  and 
Dziewonski,  1987]  and  free  oscillations  [eg.  Ritzwoller,  et  al,  1986;  Giardini,  et  al,  1987],  geoid  model¬ 
ling  [eg.  Hager,  et.  al.,  1985;  Forte  and  Peltier,  1987;  Zhang  and  Yuen,  1987],  and  earth  rotation  and 
other  geodetic  studies  [eg.  Gwinn,  et  al,  1986;  Levine,  et  al.,  1986;  Neuberg,  et  al.,  1987;  Wahr  and  de 
Vries,  1989;  Hide,  et  al.,  1989].  In  many  cases  there  is  a  possible  ambiguity  in  the  interpretation  of  the 
observations:  do  the  data  reflect  the  effects  of  core/mantle  boundary  topography  alone,  or  are  there  also 
effects  of  structure  within  the  core  or  lower  mantle?  Often,  the  possibility  of  lateral  structure  inside  the 
core  is  simply  ignored.  (A  comment  on  notation:  the  earth  is  primarily  an  ellipsoid,  with  ellipticity 
determined  by  the  hydrostatic  response  of  the  earth  to  the  centrifugal  force  associated  with  its  rotation. 
Surfaces  of  constant  density  and  boundary  surfaces  for  a  hydrostatic  earth  have  the  shape 
r  =  r0(l  +  e(r0)>,2  (©.X.)),  where  r0  is  the  mean  radius  of  the  surface,  Y$ (0,X)  is  a  spherical  harmonic, 
and  0  and  X  are  the  co-latitude  and  eastward  longitude,  respectively.  The  dimensionless  function,  e(r0), 
is  determined  by  solving  Clairaut’s  differential  equation  -  see,  for  example,  Jeffreys  [1976]. 
Throughout  this  chapter,  the  word  "hydrostatic"  will  refer  to  this  contribution  to  the  earth’s  lateral 
stratification.  The  words  "structure"  and  "topography"  will  denote  the  departure  of  a  constant  density 
surface  and  of  a  boundary,  respectively,  from  the  hydrostatic  shape.) 

Stevenson  [1987]  concludes  that  surfaces  of  constant  material  properties  inside  the  fluid  core  should 
closely  coincide  with  surfaces  of  constant  potential,  and  his  results  have  occasionally  been  invoked  to 
justify  the  assumption  of  no  internal  structure.  But,  there  is  no  guarantee  that  the  constant  potential  sur¬ 
faces  are,  themselves,  free  of  structure.  In  fact,  the  structure  of  constant  potential  surfaces  inside  the 
core  should  partly  reflect  the  topography  at  the  corc/mantle  boundary.  More  precisely,  the  constant 
potential  surfaces  inside  the  fluid  core  can  be  modelled,  as  described  below,  from  knowledge  of  the 
constant  potential  surfaces  at  the  mean  radii  of  the  core/mantle  and  inner  core/outcr  core  boundaries. 
And,  because  of  the  density  contrast  across  those  boundaries,  the  potential  at  a  boundary  radius  is 
affected  by  the  boundary  topography. 

Depending  on  the  type  of  observation  and  on  the  amplitude  of  the  internal  structure,  the  results  of 
ignoring  that  structure  when  inverting  for  core/manlle  boundary  topography  could  be  significant  For 
example,  Wahr,  et  al.  [1987]  inverted  core-phase  travel  times  to  infer  the  long  wavelength  spherical 
harmonic  coefficients  of  the  boundary  topography  for  two  cases.  In  one  case,  structure  inside  the  core 
was  ignored,  which  is  equivalent  to  assuming  that  surfaces  of  constant  potential  inside  the  core  coincide 
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with  the  hydrostatic  surfaces.  In  the  other  case,  the  core/mantle  boundary  was  assumed  to  be  a  surface 
of  constant  potential.  In  the  latter  case,  the  structure  of  the  internal  equipotential  surfaces  depends 
directly  on  the  amplitude  of  the  boundary  topography  (see  section  2  below),  and  can  be  large  for  large 
topography.  The  boundary  topography  obtained  by  Wahr,  el  al.  using  core-refracted  phases  was  substan¬ 
tially  different  for  these  two  cases  at  the  very  longest  wavelengths:  ie.  for  the  y/"(0,X)  spherical  har¬ 
monic  coefficients  for  /  <  3.  The  coefficients  were  larger,  sometimes  by  several  times,  for  the  case 
when  the  internal  structure  was  included.  The  assumption  that  the  boundary  is  an  equipotential  surface 
is  likely  to  cause  an  overestimate  of  the  internal  structure  (see  below).  Still,  the  effects  of  the  assumed 
structure  were  large  enough,  particularly  for  inversions  of  PKKP  data,  to  suggest  that  the  effects  of 
more  moderate  internal  structure  could  also  be  important. 

It  is  easy  to  understand  why  the  amplitudes  inferred  from  travel  time  studies  are  apt  to  be  larger 
when  structure  is  included.  Suppose  the  core/mantle  boundary  is  elevated  at  some  point  Q.  Then,  a 
PKP  ray  passing  through  the  boundary  at  Q  has  an  increased  travel  time  because  fast  mantle  material  in 
the  vicinity  of  Q  has  been  replaced  by  slower  core  material.  If  internal  core  structure  is  present,  sur¬ 
faces  of  constant  wave  velocity  in  the  core  beneath  Q  would  likely  be  elevated,  as  well.  The  effect  of 
this  internal  structure  is  that  core  material  at  any  level  beneath  Q  has  been  pushed  upwards,  and 
replaced  by  the  slightly  faster  core  material  immediately  below.  This  tends  to  decrease  the  travel  time 
of  the  ray,  thus  partially  offsetting  the  effects  of  the  boundary  topography.  If  the  effects  of  internal 
structure  are  not  included  in  the  inversion,  the  inversion  algorithm  interprets  the  reduced  travel  time 
perturbations  as  due  to  smaller  topography. 

The  effects  of  the  induced  internal  structure  have  been  found  to  be  less  important  for  free  oscillation 
inversions  (J.  Woodhouse,  personal  communication;  G.  Masters,  personal  communication).  The  reason 
for  this  difference  is  probably  that  the  free  oscillation  kernels  (the  functions  that  represent  the  depen¬ 
dence  of  the  eigenfrequencies  on  structure)  are  oscillatory  functions  of  radius  inside  the  core.  Since,  as 
we  shall  see,  induced  structure  inside  the  core  tends  to  have  the  same  sign  at  all  radii,  the  effects  on  the 
eigenfrequencies  would  tend  to  cancel.  On  the  other  hand,  core-phase  travel  time  kernels  have  the 
same  sign  everywhere:  an  increase  in  wave  speed  anywhere  along  the  path  reduces  the  travel  time. 
Thus,  the  effects  of  core  internal  structure  at  different  radii  tend  to  add  constructively. 

In  Section  2,  we  will  describe  the  relevant  equations,  and  their  solutions,  for  the  problem  of  comput¬ 
ing  the  structure  inside  the  core.  In  Section  3,  we  will  discuss  the  implications  of  our  results  for  the 
problem  of  inferring  core/mantle  boundary  topography  from  observations  of  the  luni-solar  nutations  and 
diurnal  earth  tides.  The  results  from  those  observations  provide  a  constraint  on  the  the  core’s  moments 
of  inertia  which,  in  turn,  can  be  related  to  the  Y”?  coefficients  of  the  core/mantle  boundary  topogra¬ 
phy  and  of  the  structure  within  the  core.  We  will  describe  the  sensitivity  of  the  results  to  different 
assumptions  about  the  internal  structure.  We  will  discuss  the  possible  effects  of  a  thin,  low  density 
fluid  layer  at  the  top  of  the  core,  such  as  might  result  from  gravitationally-induced  convection  in  the 
core.  Section  4  summarizes  the  results,  and  includes  a  discussion  of  whether  there  is  likely  to  be 
significant  internal  structure.  We  conclude  that  significant  structure  is  not  implausible,  and,  at  the  very 
least,  it  should  not  be  dismissed  out-of-hand  when  confronted  with  unexpected  seismic  anomalies. 

2  THEORY 

Stevenson  [1987]  concludes  that  dynamical  forces  originating  within  the  fluid  core  are  incapable  of 
supporting  any  appreciable  internal  structure.  Thus,  if  sizable  structure  exists,  it  must,  instead,  be  due 
to  external  forcing  from  the  mantle  or  inner  core.  In  this  section,  we  describe  a  method  for  computing 
the  structure,  given  knowledge  of  the  external  forcing. 

There  are  at  least  two  possible  ways  of  approaching  this  problem.  One  is  to  start  with  the  result  for 
a  static  fluid  that,  within  the  fluid  core,  surfaces  of  constant  density  are  surfaces  of  constant  potential. 
This  leads  directly  to  the  construction  of  an  integral-differential  equation  for  the  spherical  harmonic 
coefficients  of  the  structure  which,  after  suitable  differentiation,  reduces  to  Clairaul’s  differential  equa¬ 
tion.  This  is  the  usual  method  of  deriving  Clairaul’s  equation  [see,  for  example,  Jeffreys,  1976]. 

The  second  method  uses  the  Eulcrian  equations  of  motion  for  a  sialic  fluid  to  derive  a  differential 
equation  directly.  Boundary  conditions  are  also  obtained.  The  differential  equations  derived  using 
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these  two  methods  are  identical.  But,  with  the  first  method  it  is  the  integral-differential  equation  that 
determines  the  unknown  constants  in  the  general  solution  to  Clairaut’s  equation;  in  the  second  method  it 
is  the  boundary  conditions.  Here,  we  will  use  the  second  method. 


2a.  Forcing  From  the  Mantle  and  Core/Mantle  Boundary 

There  are  two  ways  in  which  the  mantle  can  perturb  the  constant  density  surfaces  inside  the  fluid 
core.  First,  any  non-hydrostatic  core/mantle  boundary  topography  can  cause  deformation  inside  the 
fluid  core  by  constraining  the  outer  fluid  boundary  to  have  a  given  shape.  We  will  denote  the  increase 
in  the  radial  coordinate  of  the  boundary  at  the  angular  position  (0,X)  by 

h"(QX)  =  a  £  C  >7W)  (1) 

IjH 


where  a  is  the  mean  core/mantle  boundary  radius.  Throughout  this  chapter,  the  complex  spherical  har¬ 
monics,  Y”,  are  assumed  to  be  normalized  so  that 


rr(8,X)yp'(e,X)  sine  de  dx  =  h,r  8 


m/n' 


(2) 


Second,  non-hydrostatic  mass  anomalies  in  the  mantle  can  produce  a  gravitational  force  inside  the 
core  which  deforms  the  constant  density  fluid  surfaces.  These  mass  anomalies  would  include  any 
lateral  variation  in  the  density  field  inside  the  mantle,  and  any  anomaly  due  to  a  change  in  the  shape  of 
the  core/mantle  boundary.  For  example,  if  the  boundary  is  elevated  in  some  region,  then  light  mantle 
material  in  that  region  has  been  replaced  by  heavier  core  material,  and  so  there  is  a  net  positive  surface 
mass  anomaly.  If  the  elevation  of  the  boundary  is  hM(Q,X),  then  the  effective  surface  mass  density  is: 

a(0 A)  =  (p(a >p(a+)>  hM  (0A)  (3) 

at  r=a,  where  a+  and  a_  represent  the  values  of  r  just  outside  and  just  inside,  respectively,  the 
core/mantle  boundary,  and  p(a+)  and  p(a_)  arc  the  mantle  and  fluid  core  densities,  respectively,  at  the 
boundary. 

To  describe  the  effects  of  the  mass  anomaly  inside  the  mantle  (ie.  above  the  core/mantle  boundary), 
we  define  the  gravitational  potential  energy  per  unit  mass  caused  by  the  non-hydrostatic  density  field 
inside  the  mantle  as 


V^OU)  =  £  V?m(r)  YM*) 
Is* 


(4) 


2b.  Forcing  From  the  Inner  Core  and  Inner  Core/Outer  Core  Boundary 

There  can  be  similar  forcing  of  the  fluid  outer  core  by  the  solid  inner  core.  We  describe  this  forcing 
in  terms  of  the  inner  core/outer  core  boundary  topography 

hlc{QX)  =  b  £*/$,  mW  (5) 

l,m 


where  b  is  the  mean  radius  of  the  inner  core/outer  core  boundary,  and  in  terms  of  the  gravitational 


35 


potential  energy  per  unit  mass 


V,c  (r  ,0,X)  =  I  V/Cm(r)  y,"(0A)  (6) 

IM 

caused  by  the  non-hydrostatic  density  distribution  inside  the  inner  core  (ie.,  beneath  the  inner  core/outer 
core  boundary). 


2c.  The  Response  of  the  Fluid  Core 

How  does  an  initially  spherically-symmetric  fluid  respond  to  VM ,  V,c ,  h* ,  and  h'cl  Let 


4>(r.e,*)  =  x  yr(BM 

ljn 


(7) 


be  the  total  gravitational  potential  in  the  fluid.  So, 

=  V*m  +  V<%  +  vf%  (8) 

where  the  Vf%  are  the  spherical  harmonic  coefficients  of  the  potential  caused  by  the  perturbed  density 
inside  the  fluid  core  and  by  the  surface  mass  anomalies  due  to  the  perturbed  inner  corc/outer  core  and 
core/mantle  boundaries  (for  example,  due  to  0(8, X)  in  (3)  for  the  perturbed  core/mantle  boundary). 

If  we  can  find  (<p,  „  is  unknown  because  the  fluid  core  contributions  to  V,°£  arc  unknown),  we 
can  infer  the  structure  of  constant  density  surfaces  inside  the  fluid  from  the  static-limit  requirement  that 
constant  density  surfaces  coincide  with  surfaces  of  constant  potential.  This  requirement  implies  that  if 


r 


r0  1  +  X  e^(r0)  YMX) 

Ijtl 


(9) 


describes  the  constant  density  surface  with  mean  radius  r=r0,  then 


ttjni.r)  =  - 


rg(r) 


where  g(r)  is  the  unperturbed  gravitational  acceleration. 

Dahlen  (1974;  equation  6.2]  finds  that  in  the  sialic  limit,  $/fm  is  a  solution  to: 


3r<»l .m  +  ~ 


/(/+1)  4rtG 

- 2~~  +  ~ T~  °'P 

r 2  8 


=  0 


(10) 


(11) 


where  p=p(r)  is  the  unperturbed  density.  Equation  (II)  is  the  differential  equation  for  internal  struc¬ 
ture.  It  is  equivalent  to  Poisson’s  equation  (V2$  =  4jrG8p),  where  the  source  density  (Sp)  is  the  pertur¬ 
bation  in  the  core’s  density  distribution  (use  (10)  in  (30),  below).  Clairaul’s  differential  equation  for 
e(-m(r)  be  obtained  by  using  (10)  in  (11). 

The  <pi  m  also  satisfy  internal  and  external  boundary  conditions.  At  internal  boundaries 
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(12) 


0;.*  (r~)  =  0/.*  (''+) 


4  nG  p  (r_) 
g(r_) 


4nGp(rJ 

g(r+) 


where  r_  and  r+  represent  the  values  of  r  just  inside  and  outside  the  boundary  [Dahlen,  1974;  equation 
(63)]. 


To  find  the  external  boundary  conditions 
boundary  (r  =  a)  into  the  mantle,  then  both 


note  that  if  d> i,*(r)  is  continued  across 


0/.«(O 


and 


3,0/,* (O  +  4rtGp(r)  r 


the  core/mantle 
are  continuous 


across  that  boundary  [Dahlen,  1974,  equation  (52)1.  These  continuity  conditions  can  be  replaced  by  a 
single,  inhomogeneous  boundary  condition  at  r  =  a  for  0(^,  just  below  the  boundary,  as  follows.  Note 
that  V\ f%,  Vfc,.  Vlcm,  3,  vfr, ,  and  3r  V/%  are  all  continuous  across  r=a.  (3,  V',*',  would  not  be  continu¬ 
ous  if  we  had  defined  vfin  to  include  the  effects  of  the  density  anomaly  associated  with  core/mantle 
topography.)  Also,  3,  Vf%  satisfies 


+  4«G[p(fl.)-p(a*)]  a  =  3 rVfi(a.) 


(13) 


drVf%  is  discontinuous,  because  we  have  defined  Vf%  so  that  it  includes  the  potential  from  the  per¬ 
turbed  core/manilc  and  inner  corc/outer  core  boundaries. 

By  using  and  noting  from  the  radial  dependence  of  solutions  to  Laplace's 

equation  that 


3,V£, (r)  =  jVj*. (r)  for  r  <  a 

ZMUr)  =  -  ^Vlcm{r)  for  r>b 

for  r  >  a 


(14) 


we  can  transform  (13)  to 


3,0/.*  (a)  +  —0i.« (a) : 

a 


^  V'Iw*(a)-47iG[p(a_)-p(a+) 


a  hM 


(15) 


Here,  3,0/,m(a)  is  the  value  of  3 ,0;,*(r)  just  below  the  core/mantle  boundary  (  0(im  and  are  con¬ 
tinuous  across  r-a ,  and  so  arc  the  same  on  both  sides  of  the  boundary).  A  similar  derivation  for  the 
inner  core/outer  core  boundary  gives; 


3,0/* (&)  "  Vfab)  -  4ti G 


p(fe+)-p(b_) 

>.  . 


b  hj 


ic 


06) 


where  p(ft+)  and  p(fe.)  arc  the  outer  core  and  inner  core  densities,  respectively,  at  the  inner  core/outer 
core  boundary,  and  3 r0(*(fc)  is  the  value  of  3,0/m  just  above  the  inner  core/outer  core  boundary. 

Thus,  0/^  in  the  fluid  core  is  found  by  solving  the  ordinary  differential  equation  (11),  together  with 
the  external  boundary  conditions,  (15)  and  (16),  and  using  the  continuity  conditions.  (12),  at  any  inter¬ 
nal  boundary.  Note  that  the  forcing  parameters,  v&.  v‘%>  a"*.  and  V*.  are  present  in  the  external 
boundary  conditions  but  not  in  the  differential  equation  or  the  internal  continuity  conditions. 
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The  boundary  conditions  (15)  and  (16)  suggest  that  the  fluid  core  structure  has  the  same  radial  depen¬ 
dence  for  forcing  from  the  boundary  topography  term,  hf as  from  the  applied  potential  coefficient 
V"m  (similarly  for  hlcm  and  !'/$.)•  The  reason  is  that  the  forcing  from  h^,,  like  the  forcing  from  Vfm, 
is  actually  gravitational.  By  fixing  the  perturbed  surface  of  the  core  to  have  the  shape  specified  by  hM , 
we  are  introducing  the  effective  surface  mass  density,  (3),  at  r=a,  which  acts  gravitationally  on  the 
material  inside  the  fluid  core  and  deforms  it. 

In  some  cases  there  may  be  more  useful  ways  to  parameterize  the  forcing  than  in  terms  of  the  boun¬ 
dary  topography  and  applied  gravitational  potential  coefficients.  One  particularly  simple  alternative  is 
to  use,  instead,  the  Y”  coefficients  of  the  equipotentiai  surface  at  the  inner  and  outer  core  boundaries  as 
the  forcing  parameters.  A  problem  with  this  alternative  is  that  the  shape  of  an  equipotentiai  surface 
depends  on  the  structure  inside  the  fluid  core,  which  is  assumed  here  to  be  unknown.  Still,  as  we  shall 
see,  this  choice  of  parameterization  can  be  convenient  for  physical  interpretation  and  is  particularly  use¬ 
ful  when  inverting  for  the  structure. 

To  include  this  parameterization,  let 


=  a^eLm(a)Yr(QX) 

Ijn 

=  6Le,,„(f>)>7W) 

l  Si 


(17) 


represent  the  increase  in  the  radial  coordinates  of  the  equipotentiai  (and  equidensity)  surfaces  at  the 
core/mantle  and  inner  core/outcr  core  boundaries,  respectively.  Then,  the  boundary  conditions  (15)  and 
(16)  are  replaced  by  the  much  simpler  conditions  (using  (10)) 


4>/,n(tf)  =  -  ag(a)e,_m(a)  (18) 

*'„(!>)  =  -bk(b)zl/H(b)  (19) 


2d.  Approximate  Solutions 

Suppose  the  fluid  core  was  homogeneous,  with  uniform  density  p(r)  =  pQ.  Then,  the  general  solution 
to  the  differential  equation  (11)  would  be 


Qljn 


<!  +  !> 


(20) 


so  that,  using  (10)  with  g(r)  =  y7tGp</  (we  assume 


the  inner  core  also  has  uniform  density  p0), 


r  i 

1-2 

'  , 

E/, n(r)  =  F 

k 

i  a 

+  G 

a 

r 

(21) 


where  the  constants  F  and  G  are  proportional  to  D  and  E . 

The  solution  for  an  inhomogeneous  fluid  core  is,  of  course,  more  complicated.  We  have  solved  (11) 
numerically  using  p(r)  from  earth  model  PREM  [Dzicwonski  and  Anderson,  1981),  and  then  using  per¬ 
turbations  of  PREM  obtained  by  inserting  a  thin,  low  density  layer  at  the  top  of  the  fluid  core,  such  as 
might  result  from  compositional  convection  in  the  core  (see  section  3,  below,  for  the  motivation  for 
including  such  a  boundary  layer).  We  find,  in  all  these  cases,  that  e/^,(r)  has  very  close  to  the  same 
radial  dependence  as  docs  the  homogeneous  solution  (21).  The  gravitational  potential,  ,  agrees  less 
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well  with  the  homogeneous  result  (20).  But,  by  using  (21)  in  (10),  we  can  obtain  a  reasonably  accurate 
approximation  for  <}>,  m: 


1 

'  ' 

1-2 

f  1 

(l*2T 

1 

II 

5 

F 

r 

a 

+  G 

7 

l/j 

g{r)r 


(22) 


where  F  and  G  are  constants  that  must  be  determined  from  the  boundary  conditions.  Using  (22)  either 
in  (15)  and  (16)  or  in  (18)  and  (19),  depending  on  the  type  of  boundary  conditions  being  considered, 
and  solving  for  F  and  G ,  we  obtain  results  for  the  internal  structure  coefficient  £/im(r)  that  agree  well 
with  the  more  accurate  numerical  solution  to  (11).  For  example,  we  find  that  for  boundary  conditions 
(18)  and  (19)  with  Zi^,(a)  =  1  and  £/.„(£)  =  0  (imposed  structure  at  the  core/mantle  boundary  and  none 
at  the  inner  core/outcr  core  boundary),  and  for  every  /  <  10,  the  resuit  for  £/m  (r)  at  any  radius,  r, 
derived  using  the  approximate  form  (21)  agrees  with  the  numerical  solution  using  PREM  to  better  than 
2%  of  the  value  of  £/-m(a)  (e,  „(a)  is  the  maximum  value  of  ttjK(r)  for  these  boundary  conditions). 
The  agreement  for  (18)  and  (19)  with  e/m(a)  =  0  and  e(-m(h)  =  1  (imposed  structure  at  the  inner 
core/outer  core  boundary)  is  better  than  11%  of  £/„(&)  (£(|M(6)  is  the  maximum  value  of  zt/n(r)  in 
this  case). 

As  a  special  case  which  will  be  considered  in  more  detail  in  the  next  section,  suppose  the  inner  core 
behaves  like  a  fluid  over  long  time  periods.  Then,  we  can  extend  our  model  and  compute  the  response 
of  the  entire  (outer  plus  inner)  core  to  forcing  from  the  mantle  alone.  In  that  case,  <j>,  m  must  satisfy  the 
differential  equation  (11)  for  0  S  r  <  a,  it  must  be  finite  when  r=0,  and  it  needs  to  satisfy  boundary 
conditions  only  at  r-a :  cither  (15)  or  (18).  For  either  boundary  condition,  the  homogeneous  solutions 


=  F 


0  ,„(r)  =  -F 


r  1 1-2 
r_ 

a 

i-i 


(23) 


g(r)a 


are  found  to  closely  approximate  the  solutions  to  the  differential  equation  (11)  for  realistic  p(r).  For 
boundary  condition  (15), 


F  = 


4wGp(fl_)  a  h &  ~  Vt!m(o){2l+l)la 
2(l-\)s(a)  +  4nGap(a.) 


(24) 


For  boundary  condition  (18): 


(25) 


These  homogeneous  results  are  identical  to  the  exact  solutions  for  1= 1  no  matter  what  the  functional 
form  of  p(r).  For  1  <  /  <  10,  they  provide  estimates  for  e(m  (r)  which  are  accurate  to  better  than  2% 
of  tijnia)  (e(  m(a)  is  the  maximum  value  of  £/^(r)  in  this  case). 

The  result  (24)  for  F  implies  that  the  response  of  the  core  to  (a )  increases  with  decreasing  p(a_). 
As  an  extreme  case,  £,.,^n(r)  is  inversely  proportional  to  p(a_).  This  behavior  is  readily  understood 
physically,  as  follows.  A  V"]jn(a)  term  represents  a  uniform  force  on  the  core.  Because  of  this  force, 
the  surfaces  of  constant  density  inside  the  core  are  all  displaced  the  same  amount  along  the  direction  of 
the  force,  as  shown  in  Figure  IV.  1.  This  displacement  causes  the  layer  at  the  outer  surface  of  the  core 
to  develop  the  asymmetrical  shape  shown  in  Figure  IV.  1:  extra  mass  on  one  side,  a  reduction  in  mass 
on  the  other.  This  mass  asymmetry  provides  an  additional  gravitational  force  on  the  core  interior. 
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which  opposes  the  force  from  V{iljn.  The  displacements  in  the  core  increase  until  these  two  forces 
exactly  cancel.  Since  the  asymmetry  in  the  surface  mass  distribution  is  proportional  to  p(aJef,tiM(<i), 
and  since  the  gravitational  force  from  the  layer  is  proportional  to  the  mass  asymmetry,  then  e(=1-m(a)  is 
inversely  proportional  to  p(a_).  The  dependence  of  ztjH (a)  on  p(a_)  for  />  I  can  be  understood  simi¬ 
larly. 


3  THE  NUTATION  AND  EARTH  TIDE  CONSTRAINT 


3a.  Constraints  on  the  Core’s  Inertia  Tensor 

The  free  core  nutation  (FCN)  is  a  normal  mode  of  the  earth,  with  an  approximately  diurnal  period 
and  with  an  eigenfunction  consisting  primarily  of  a  relative  rotation  between  the  core  and  mantle.  The 
mode  has  been  modelled  semi-analytically  for  an  earth  with  slightly  aspherical,  but  otherwise  arbitrary 
structure,  and  for  an  inner  core  which  is  assumed  to  be  fluid  [see,  eg.  Wahr  and  dc  Vries,  1989].  The 
FCN  eigenfrequency  in  that  case  is  found  to  depend  on  the  dynamical  ellipticity  of  the  core: 

e  =  (C -A  )IA  (26) 

where  C  and  A  are  principal  moments  of  inertia  of  the  core.  The  effects  on  the  eigenfrequency  of  a 
solid  inner  core  are  small  [de  Vries  and  Wahr,  1989],  and  need  not  concern  us,  here. 

If 


8p(r,0,X)  =  £8p  /W/W) 

Ijn 


(27) 


is  the  aspherical  density  distribution  inside  the  core,  then  the  contribution  to  e  from  the  non-hydrostatic 
structure  and  topography  is 


8e  =  - 


2VS/5 


a5/i"oP(o-)  +  |6p2V)r4dr 


(28) 


so  that  Se  depends  only  on  the  X®  components  of  boundary  topography  and  internal  structure. 

The  luni-solar  forced  nutations  and  diurnal  earth  tides  occur  at  frequencies  close  enough  to  the  FCN 
eigenfrequency  to  be  significantly  affected  by  the  mode.  Their  amplitudes  have  been  used  to  constrain 
the  eigenfrequency  (see,  eg.,  Gwinn,  et  al.,  1986,  for  nutation  results;  Levine,  et  al,  1986,  and  Ncubcrg, 
ct  al.,  1987,  for  tidal  results].  The  nutation  and  earth  tide  results  for  the  eigenfrequency  arc  consistent 
with  each  other,  and  imply  a  value  for  e  about  5%  in  excess  of  the  hydrostatic  value.  The  VLBI  nuta¬ 
tion  estimates,  which  currently  provide  the  most  accurate  results  for  the  eigenfrequency,  imply  that  the 
contribution  to  e  from  non-hydrostatic  internal  core  structure  and  core/mantle  boundary  topography  is 

Se  =  1.22  x  10-4  (29) 

with  an  uncertainty  of  about  ±  10%  [Gwinn,  et  al.,  1986].  Since 

5 pf (r  )  =  ~rt,_m{r  )dr  p(r )  (30) 


then  (28)  and  (29)  imply  the  following  constraint  on  h^o  and  z^r ): 


40 


2n'h/5 

A 


a 


a^hi.oP (a-)  -  \r5t2,0(r)drp(,r)dr 
0 


-1.22x1  O'4 


(31) 


3b.  Implications  for  Core/Mantle  lioundary  Topography 

Consider  the  situation  described  at  the  end  of  Section  2,  where  the  entire  core  responds  as  a  fluid  to 
forcing  from  the  mantle  and  corc/mantle  boundary.  Parameterizing  the  internal  structure  in  terms  of  the 
cquipotcntial  surface  components,  e,  n(a),  we  find  that 


^  (r )  =  f  2.m  (a )  =  constant  (32) 

is  a  good  approximation  to  the  1=2  numerical  solution  for  realistic  density  p(r),  particularly  for  large  r. 
The  large  r  results  arc  especially  important  here  Nscausc  the  integrand  in  (31)  includes  a  factor  of  r5, 
which  results  in  heavy  weighting  of  large  r  values.  The  r5  factor  also  permits  us  to  ignore  the  possible 
effects  of  independent  structure  in  the  inner  core. 

Using  (32)  in  the  constraint  (31),  and  integrating  by  parts  gives: 

3  , —  2  Vrt/5  p (a-)a5  C  u  3  , 

j  V 5/rr  e2.0(n )  + - j  o-e2.o(a  )J  =  -1 .22x10-*  (33) 

When  evaluated  numerically,  this  result,  (33),  gives  a  combined  constraint  on  e^a)  and  h%  which 
agrees  to  better  than  1%  with  the  constraint  obtained  using  more  accurate  numerical  solutions  for  e20(r) 
in  (31).  The  advantage  of  (33)  is  that  its  dependence  on  p(a_)  is  clear. 

Results  for  PRKM 

Using  valu^^if  p (a.)  and  A  that  arc  consistent  with  the  PREM  density  distribution,  (33)  reduces  to 
the  constraint: 


h%  +  .077  e20(a)  =  -1.39xlO-4 


(34) 


Note  that  (34)  is  a  constraint  primarily  on  the  boundary  topography  coefficient,  h20 .  There  is  little  sen¬ 
sitivity  to  the  structure  of  the  cquipotcntial  surface,  t)/n(a)  (ic.  .077  is  much  less  than  1),  because:  (1) 
the  components  of  the  core’s  inertia  tensor  (in  this  case,  C  and  A )  arc  much  more  sensitive  to  the  den¬ 
sity  distribution  in  the  outer  regions  of  the  core  than  to  the  density  distribution  in  the  interior  (the  r5 
dependence  in  (31)  );  and  (2)  PREM  has  lildc  radial  density  contrast  through  the  outer  region  of  the 
core,  and  so  it  doesn’t  much  matter  what  the  shapes  of  the  constant  density  surfaces  are.  The  con¬ 
straint  is  strongly  sensitive  to  the  topography  coefficient,  h2,o<  because  there  is  a  large  contrast  in  the 
core  density  across  the  core  boundary:  from  about  9.9  gm/cm3  just  below  the  boundary  to  0  gm/cm3 
above  (  0  gm/  cm3  because  there  is  no  core  material  above  the  core/mantlc  boundary). 

Unless  the  cquipotcntial  surface  at  r=a  has  Y2  structure  which  is  many  times  the  size  of  the  Y2 
boundary  topography,  the  result  (34)  gives  approximately 

a  hi o  =  j  km  (35) 
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Effects  of  a  thin,  core  boundary  layer 

The  constraint  on  the  boundary  topography  is  different  if  there  is  a  thin,  low  density  fluid  layer  at  the 
top  of  the  core,  such  as  might  be  caused  by  compositional  convection.  In  that  case,  the  sensitivity  of 
the  core’s  inertia  tensor  to  the  equipotential  coefficient,  e2  0(a ),  is  increased  because  of  the  larger  radial 
stratification  in  the  outer  region  of  the  core.  At  the  same  lime,  the  sensitivity  to  boundary  topography 
is  decreased  because  of  the  smaller  density  contrast  across  the  core/mantle  boundary.  In  effect,  some  of 
the  9.9  gm/cm3  core  density  contrast  that  occurs  across  the  core  boundary  in  PREM  is  transferred  to 
one  or  more  surfaces  (or  across  a  finite  region)  within  the  boundary  layer.  The  larger  the  radius  of  one 
of  those  internal  surfaces,  the  larger  the  effect  on  the  inertia  tensor  of  the  density  contrast  across  that 
surface.  Thus,  the  importance  of  the  e20(a )  term  in  the  constraint  must  depend  partly  on  the  boundary 
layer  thickness,  as  well. 

In  anticipation  of  these  results,  we  define  the  dimensionless  parameter: 


8rc  P(a-)a5 
A  =  l5~~A 


(36) 


The  core’s  principal  moment,  A ,  depends  on  the  density  and  thickness  of  the  boundary  layer.  But,  that 
dependence  is  minimal.  For  example,  suppose  the  boundary  layer  has  thickness  D.  and  suppose  the 
density  throughout  the  layer  is  smaller  than  the  PREM  density  distribution  by  an  amount  8p  =  constant. 
Then,  for  D  c  a: 


A  =  A0-Da45p87t/3 


(37) 


where  A  0  is  the  PREM  value  of  A  and  is  approximately  9.12xl043gm  cm2.  To  estimate  the  possible 
effects  of  D ,  note  that  the  density  in  the  boundary  layer  is  unlikely  to  be  less  than  about  5.6  gm/cm3, 
which  is  the  PREM  density  at  the  bottom  of  the  mantle,  and  so  5p  S  9.9  -  5.6  =  4.3  gm/cm3.  We  find 
numerically: 


A  =  A0 

j  Da45p87t/3' 

<A0 

1-2.0  — 

a 

(38) 


This  result  shows  that  as  long  as  the  boundary  layer  is  much  thinner  than  the  core  radius  (Dca), 
A  =A0,  and  so  A  is  relatively  insensitive  both  to  the  boundary  layer  density  and  to  the  thickness. 
Thus,  the  only  important  dependence  of  A  on  the  boundary  layer  is  through  p(«_)  in  the  numerator  of 
(36).  If  we  assume  that  p(a_)  is  no  larger  than  the  PREM  density  at  the  top  of  the  core,  and  no  smaller 
than  the  PREM  density  at  the  bottom  of  the  mantle,  then 


0.522  <  A  <  0.928 


(39) 


With  the  parameter  A,  (33)  reduces  to 

h%0  A  +  £2,0 (a  )(1-A)  =  -1.29x10"  (40) 


(Using  the  PREM  result,  A  =  0.928,  in  (40),  wc  obtain  (34).) 

Consider  the  following  two  examples.  First,  consider  the  extreme  case  that  the  corc/mantlc  boundary 
is  an  equipotential  surface.  Then  e,  m(a  )=/i/Mm,  and  (40)  reduces  to 
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(41) 


h%0  =  -1-29x1  O'4 

There  is  no  dependence  on  A  in  (41).  The  reason  is  that  in  this  example  the  internal  boundary  layer 
structure  has  the  same  (9,X)  dependence  as  the  core/mantlc  topography.  So,  it  doesn’t  matter  whether 
the  9.9  gm/cm3  density  contrast  is  spread  out  across  the  top  layers  of  the  core  or  is  all  concentrated 
right  at  the  boundary;  the  effects  on  the  inertia  tensor  arc  the  same. 

Second,  suppose  there  is  no  1=2,  m  =  0  structure  inside  the  core,  so  that  £2,0(0)  =  0.  Then,  (40) 
reduces  to: 

/t£0A  =  -1.29x10^  (42) 

It  is  clear  from  (42)  that  the  nutation  and  earth  tide  constraint  is  able  to  accommodate  larger  h^Q 
values  if  A  is  smaller  than  the  PREM  value.  Smaller  A  (ie.  smaller  p(a_))  implies  that  there  is  a 
smaller  density  contrast  across  the  core/mantlc  boundary,  and  so  the  boundary  topography  must  be 
larger  to  cause  a  given  perturbation  in  the  inertia  tensor.  If  A  is  equal  to  the  minimum  value  of  0.522 
(see  equation  (39)),  then  (42)  implies  that  /i^0  =  2.5x10^.  This  numerical  result  is  probably  the  largest 
possible  value  of  h -/0  that  is  consistent  with  the  assumption  of  no  structure  inside  the  core.  It 
corresponds  to  Y°  boundary  topography  of  .87  km. 

But,  this  last  example  (£2,0(0  )=0)  is  not  the  most  extreme  possibility.  The  nutation  constraint  for  a 
core  with  a  boundary  layer  can  accommodate  still  larger  boundary  topography,  if  the  internal  structural 
parameter,  £2.0(0),  is  the  same  order  as  h^o  but  with  the  opposite  sign.  For  example,  if  we  use  the 
lower  bound  A  =  0.522  in  (40),  we  obtain 


A&+.916  e2.o(0)  =  -2.47x  UT* 


(43) 


For  A  =  0.522,  roughly  half  the  9.9  gm/cm3  density  contrast  occurs  inside  the  core,  where  the  contribu¬ 
tions  to  the  inertia  tensor  depend  on  e2,o(0)  instead  of  on  h"0.  This  is  why  the  constraint  equation, 
(43),  is  as  sensitive  to  the  internal  structure  as  to  the  boundary  topography. 

Suppose  h" 0  and  £2,0(0)  have  opposite  signs.  Then,  the  effects  of  the  boundary  topography  on  the 
inertia  tensor  arc  partially  offset  by  the  effects  of  the  structure  in  the  boundary  layer,  and  both  /i^0  and 
£2,0(0 )  can  have  large  amplitudes  and  still  satisfy  the  nutation/carth  tide  constraint.  There  can  be  large 
1=2  topography  for  A  greater  ihan  0.522,  as  well.  However,  note  from  (40)  that  the  closer  A  is  to  1,  the 
larger  £2,0(0)  must  be  to  offset  the  effects  of  the  boundary  topography. 

Although  this  combination  of  a  low  density  boundary  layer  and  roughly  equal  but  opposite  internal 
structure  and  topography  could  reconcile  the  nutation  constraint  with  large  topography,  we  arc  not 
presently  advocating  such  a  situation.  For  one  thing,  there  is  currently  no  compelling  reason  to  believe 
that  the  1=2,  m=  0  boundary  topography  is  much  in  excess  of  1/2  km.  For  another,  there  are  dynamical 
considerations  that  tend  to  argue  against  such  a  configuration.  These  arguments  are  presented  at  the 
end  of  the  following  section. 

4  DISCUSSION 

There  arc  two  parts  to  this  chapter.  In  the  first  part,  we  show  that  there  can  be  non-hydrostatic, 
laterally  varying  structure  inside  the  outer  core.  Such  structure  could  be  caused  by  gravitational  forcing 
from  non-hydrostatic  mass  anomalies  inside  the  mantle  or  inner  core  or  from  the  density  anomalies 
associated  with  non-hydrostalic  core/mantlc  or  inner  corc/outer  core  boundaries. 

The  differential  equation  for  the  core  structural  parameters,  £(^,(r),  is  given  by  (10)  and  (11).  We 
consider  two  types  of  possible  boundary  conditions  for  this  problem.  The  boundary  conditions  (15)  and 
(16)  are  appropriate  if  the  core  boundary  topography  and  the  density  distribution  inside  the  mantle  and 
inner  core  are  known.  The  boundary  conditions  (18)  and  (19)  arc  appropriate  if  the  equipotential  sur¬ 
faces  at  the  mean  inner  core  and  outer  core  radii  arc  known. 
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It  is  unlikely  that  the  cquipotcntial  surfaces  at  the  boundaries  would  be  known.  After  ail,  the  poten¬ 
tial  depends  partly  on  the  density  distribution  inside  the  core  and  that  is  not  known  prior  to  solving  the 
equations.  If  the  inner  and  outer  core  boundaries  could  be  assumed  to  be  equipotential  surfaces,  then 
the  e/jn(a)  and  tlfK{b)  in  (18)  and  (19)  could  be  determined  from  knowledge  of  the  boundary  topogra¬ 
phy.  However,  this  assumption  is  unlikely  to  be  valid,  because  it  would  probably  require  unacceptably 
small  values  for  the  viscosity  of  the  lower-most  mantle  (Brad  Hager,  personal  communication). 

Instead,  the  primary  use  of  the  boundary  conditions  (18)  and  (19),  is  in  providing  a  convenient 
parameterization  for  the  internal  structure.  One  of  the  possible  applications  of  the  theory  described  here 
is  to  facilitate  data  inversions  for  the  structure.  Seismic  travel  lime  studies,  for  example,  might  be  par¬ 
ticularly  useful.  The  results  of  this  chapter  show  that  inversions  for  core  structure,  unlike  inversions  for 
mantle  structure,  do  not  require  solutions  for  the  structure  at  different  radii.  Instead,  the  radial  depen¬ 
dence  of  the  structure  inside  the  fluid  core  is  determined  by  the  differential  equation,  (11),  and  so  the 
entire  solution  can  be  parameterized  in  terms  of  a  set  of  boundary  coefficients:  the  right  hand  sides  of 
either  (15)  and  (16)  or  of  (18)  and  (19).  The  data  can  be  inverted  to  find  those  coefficients.  The  boun¬ 
dary  conditions  (18)  and  (19)  arc  particularly  useful  because  of  their  simplicity. 

The  numerical  solutions  to  the  differential  equation  (11)  arc  found  to  have  an  easily  characterized, 
approximate  radial  dependence,  given  by  (21).  For  short  horizontal  wavelengths  (large  /),  the  structure 
falls  off  rapidly  with  distance  from  the  boundaries.  For  example,  the  structural  amplitudes  beneath  the 
core/mantle  boundary  decrease  with  decreasing  radius  approximately  as  (r/a)(,~1)  (the  amplitude  of  the 
structure  is  rel/n(r),  and  el/B(r)  is  proportional  to  r<,_2):  see  equation  (23)).  This  implies  that  there  is 
not  apt  to  be  much  short  wavelength  core  structure  away  from  boundaries.  It  helps  explain  the  results 
of  Wahr,  et  al.  [1987],  who  found  that  their  seismic  travel  lime  inversions  for  corc/mantle  boundary 
topography  were  relatively  insensitive  to  assumptions  about  the  internal  structure  at  angular  orders  /  > 
3.  Note,  also,  that  {r  la )(,~"  is  a  monotonic  function  of  radius.  As  described  in  the  Introduction,  this 
helps  explain  why  fluid  core  structure  has  a  more  important  effect  on  travel  times  than  on  free  oscilla¬ 
tions. 

In  the  second  part  of  this  chapter,  we  discuss  the  possible  effects  of  fluid  core  structure  on  the 
interpretation  of  nutation  and  diurnal  earth  tide  observations.  Those  observations  constrain  the  core’s 
dynamical  ellipticity,  which  is  related  to  the  core's  principal  moments  of  inertia. 

The  moments  of  inertia  can  be  affected  by  boundary  topography  and  by  internal  structure.  The  sensi¬ 
tivity  to  the  topography  and  structure  depends  on  the  radial  density  distribution  at  the  top  of  the  core. 
If  the  density  is  consistent  with  PREM,  there  is  little  sensitivity  to  internal  structure.  The  PREM  den¬ 
sity  varies  slowly  with  radius  in  the  outer  region  of  the  core,  and  so  the  shapes  of  the  constant  density 
surfaces  are  pretty  much  irrelevant.  In  this  case,  the  results  for  the  dynamical  ellipticity  mostly  con¬ 
strain  the  core/mantle  boundary  topography. 

If  there  is  a  low  density  boundary  layer  at  the  top  of  the  core,  the  inertia  tensor  is  more  sensitive  to 
internal  structure  and  less  sensitive  to  the  boundary  topography.  To  understand  this  result,  note  that  the 
outer  surface  of  the  boundary  layer  would  coincide  with  the  core/mantle  boundary,  while  the  shape  of 
the  boundary  layer  inner  surface  would  depend  on  the  internal  structure.  The  smaller  the  density  of  the 
boundary  layer,  the  larger  the  density  contrast  across  the  inner  surface,  and  so  the  greater  the  sensitivity 
to  internal  structure.  At  the  same  lime,  a  smaller  boundary  layer  density  means  decreased  sensitivity  to 
the  shape  of  the  outer  surface  of  the  boundary  layer:  ic.  to  the  core/mantle  boundary  topography. 

The  general  result  for  the  nulation/earth  tide  constraint  is  given  by  (40).  The  dimensionless  parame¬ 
ter  A  in  that  equation  depends  on  the  density  at  the  core/mantle  boundary  (sec  equation  (36)),  and  is 
unlikely  to  have  a  numerical  value  outside  the  interval  given  in  (39).  Small  values  of  the  boundary 
layer  density  imply  small  values  of  A.  The  PREM  density  distribution  (A  =  0.928),  leads  to  an  inferred 
y®  coefficient  of  the  core/mantle  boundary  topography  of  about  1/2  km. 

A  larger  numerical  result  for  the  topography  is  obtained  if  there  is  a  low  density  boundary  layer  (A  < 
0.928)  at  the  top  of  the  core.  Still,  using  the  lower  bound  for  A  of  0.528,  we  find  that  the  only  way  the 
K®  topography  coefficient  could  be  larger  than  .87  km,  is  if  there  is  internal  structure  with  the  opposite 
sign  of  the  boundary  topography.  In  that  case,  contributions  to  the  moments  of  inertia  from  the  topog¬ 
raphy  are  partially  offset  by  contributions  from  the  structure  at  the  bottom  of  the  boundary  layer. 
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As  an  example  of  how  these  constraints  can  be  used,  suppose  we  invert  a  seismic  data  set  to  leant 
about  core/mantle  boundary  topography,  and  obtain  an  1=2,  m= 0  coefficient  with  an  amplitude 
significantly  larger  than  1/2  km.  Our  result,  then,  is  at  odds  with  the  nutation  and  earth  tide  constraint 
estimated  using  the  PREM  value  of  A.  To  reconcile  our  seismic  coefficient  with  that  constraint,  we 
might  be  tempted  to  hypothesize  the  existence  of  a  low  density  boundary  layer  at  the  top  of  the  core, 
and  of  internal  structure  with  a  large  amplitude  but  with  the  opposite  sign  of  the  boundary  coefficients. 

One  problem  we  might  encounter,  though,  is  that  the  internal  structure  we  arc  proposing  could  be 
inconsistent  with  the  assumptions  we  made  before  we  started  our  seismic  inversion.  And,  the  effects  on 
our  inversion  of  ignoring  this  structure  could  be  important.  For  example,  the  1=2  topography 
coefficients  inferred  by  Wahr,  et  al.  [1987]  from  PKP  travel  times,  differed  by  up  to  a  factor  of  two 
depending  on  whether  they  assumed  the  core/mantle  boundary  was  an  equipotential  surface  or  ignored 
internal  structure  entirely. 

In  fact,  even  the  hypothesis  of  a  boundary  layer  itself,  without  specifying  any  internal  structure,  could 
have  consequences  for  our  inversion.  For  example,  the  sensitivity  of  PKP  travel  times  to  corc/mantlc 
boundary  topography  depends  on  the  velocity  contrast  across  the  boundary.  And  a  core-side  boundary 
layer  could  affect  that  contrast. 


Is  There  Likely  to  be  Significant  Internal  Structure? 

An  even  more  fundamental  issue,  is  whether  the  combination  of  large  opposing  structure  and  topogra¬ 
phy  and  a  low  density  fluid  boundary  layer  is  plausible,  based  on  our  understanding  of  mantle  dynam¬ 
ics.  In  fact,  the  question  of  whether  there  could  be  any  appreciable  structure  inside  the  core  is  clearly 
important  for  all  applications  of  the  results  described  in  this  chapter,  whether  the  core  has  a  low  density 
boundary  layer,  or  not. 

What  insight  we  can  offer  into  this  problem  is  based  on  preliminary  results  from  Dehant  and  Wahr 
(in  preparation),  who  compute  the  deformation  of  internal  boundaries  and  the  perturbation  in  the  gravi¬ 
tational  potential  caused  by  internal  loading,  using  methods  similar  to  those  employed  by  Richards  and 
Hager  [1984],  Forte  and  Peltier  [1987],  and  Zhang  and  Yuen,  [1987],  The  results  are  strongly  depen¬ 
dent  on  the  assumed  characteristics  of  the  D”  layer  at  the  base  of  the  mantle,  and  on  whether  there  is  a 
low  density,  fluid  boundary  layer  al  the  top  of  the  core.  For  example,  if  convection  in  the  lower  mantle 
penetrates  all  the  way  to  the  base  of  the  mantle  (ic.  to  the  corc/mantlc  boundary),  if  the  viscosity  of  the 
D”  layer  is  not  very  much  smaller  than  the  viscosity  of  the  rest  of  the  lower  mantle,  and  if  there  is  no 
core-side  boundary  layer,  then  the  /  =2  structural  coefficients  (e2,m  (a ))  are  roughly  20%  of  the  /  =2  topo¬ 
graphic  coefficients  (/t^),  for  loads  in  the  mid-mantle. 

On  the  other  hand,  if  there  is  a  barrier  to  convection  somewhere  in  the  D”  layer  above  the 
corc/mantle  boundary,  if  the  D”  layer  has  a  significantly  lower  viscosity  than  the  rest  of  the  lower  man¬ 
tle,  and  if  there  is  no  low  density  fluid  boundary  layer  at  the  top  of  the  core,  then  a  mid-mantle  load 
could  cause  perturbations  in  £2,m(a)  and  h%m  of  roughly  equal  amplitudes  but  opposite  signs.  (The 
opposing  signs  of  the  boundary  and  equipotential  surface  in  this  case  are  analogous  to  a  similar  effect  at 
the  outer  surface,  discussed  by  Richards  and  Hager  [1984],  that  occurs  when  there  is  a  barrier  to  flow 
between  the  upper  and  lower  mantles.)  We  note  that  the  case  of  roughly  equal  amplitudes  occurs 
because,  for  an  internal  load  of  given  strength,  the  topography  becomes  small,  and  not  because  the 
internal  structure  becomes  large. 

Thus,  it  is  possible  that  the  induced  structure  could  be  as  large  or  larger  than  the  induced  topography. 
This  would  occur,  for  mid-mantle  loads,  for  a  restricted  but  not  implausible  range  of  D”  parameter 
values.  This  suggests  that  the  possibility  of  significant  internal  structure  should  not  be  dismissed  out- 
of-hand  when  interpreting  unexpected  seismic  anomalies. 

If  a  low  density  boundary  layer  is  included  at  the  top  of  the  core,  the  induced  boundary  topography  is 
increased.  The  boundary  tends  to  deform  in  response  to  an  internal  load,  so  that  the  total  surface  mass 
anomaly  associated  with  the  deformation  equals  some  value  that  is  determined  by  the  amplitude  and 
wavelength  of  the  load.  The  surface  mass  anomaly,  (3),  is  the  product  of  the  deformation  ( hM )  and  the 
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density  contrast  across  the  boundary  (p(a.)-p(a+)).  If  there  is  a  low  density  boundary  layer  at  the  top 
of  the  core,  then  the  density  contrast  is  decreased,  and  so  the  topography  must  be  correspondingly 
increased. 

This  argues  against  the  combination  of  a  low  density  boundary  layer  and  roughly  equal  but  opposite 
core  structure  and  topography,  introduced  at  the  end  of  Section  3  as  a  way  of  reconciling  large  boun¬ 
dary  topography  with  the  nutation  constraint.  If  the  boundary  layer  density  is  low  enough  that  there  is 
only  a  small  density  contrast  across  the  core/mantle  boundary,  then  the  load-induced  topography  on  the 
boundary  would  almost  certainly  be  much  larger  than  the  load-induced  internal  structure.  We  thus  con¬ 
clude  that  the  nutation  constraint  is  probably  not  compatible  with  1=  2,  /n=0  boundary  topography  much 
in  excess  of  .87  km  (the  result  for  A  =  .528  and  no  internal  structure).  Furthermore,  we  emphasize  that 
there  is  currently  no  compelling  reason  to  believe  that  the  1=2,  m  = 0  topography  is  anything  different 
than  the  1/2  km  result  implied  by  the  nutation/body  tide  observations  using  the  PREM  value  of  A  = 
.928. 
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CHAPTER  IV  FIGURE  CAPTIONS 


Figure  IV.  1.  An  illustration  of  the  response  of  the  cote  to  an  applied  /  a  1  potential.  The  potential 
induces  a  uniform  translation  of  the  constant  density  surfaces  inside  the  core.  This  causes  the  layer 
at  the  top  of  the  core  to  have  an  aspherical  shape,  with  extra  mass  on  one  side  and  a  reduction  of 
mass  on  the  other.  This  density  anomaly  causes  a  gravitational  force  on  the  core  interior  which 
opposes  the  applied  force. 
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CHAPTER  V 


A  COMPARISON  OF  NMC  AND  FNOC  SEA  LEVEL  PRESSURE  VALUES 


SUMMARY 


A  comparison  between  values  of  sea  level  atmospheric  pressure  provided  by  the  National  Meteoro¬ 
logical  Center  and  the  Fleet  Numerical  Oceanographic  Center  was  performed  for  time  periods  coin¬ 
cident  with  the  SEASAT  and  GEOSAT  Altimeter  Missions.  Global  sea  level  pressure  grids,  each 
corresponding  to  one  of  three  3-day  repeat  orbits  of  the  satellite,  were  generated  for  a  spatial  and  tem¬ 
poral  comparison.  The  results  indicate  that  differences  greater  that  10  mbar  are  quite  common  in  the 
southern  ocean  with  offsets  as  large  as  40  mbar  in  the  Southern  Pacific. 


1.  INTRODUCTION 


To  interpret  satellite  altimeter  data  in  terms  of  wind  generated  currents,  it  is  first  desirable  to  remove 
pressure-driven  variations  in  sea  level  from  the  data.  This  is  usually  achieved  by  assuming  an  inverted 
barometer  response,  which  postulates  an  inverse  relationship  between  sea  level  atmospheric  pressure  and 
sea  surface  height  (-l.Olcm/mbar).  To  apply  this  correction,  it  is  necessary  to  adopt  some  model  for 
sea  level  pressure  over  the  global  ocean  as  a  function  of  time.  For  the  SEASAT  and  GEOSAT  altime¬ 
ter  data,  along  track  estimates  of  the  pressure  are  obtained  by  interpolating  between  global  sea  level 
pressure  grids  provided  by  the  Fleet  Numerical  Oceanographic  Center  (FNOC). 

The  accuracy  of  these  pressure  corrections  depends  critically  on  the  accuracy  of  the  FNOC  pressure 
data.  In  this  chapter  we  compare  those  data  with  pressure  fields  provided  by  the  National  Meteorologi¬ 
cal  Center  (NMC)  for  a  time  periods  corresponding  to  the  SEASAT  and  GEOSAT  Altimeter  Missions. 
Spatial  and  temporal  differences  between  the  data  sets  are  calculated  and  plotted. 


2.  METHOD 


The  FNOC  data  used  in  this  comparison  were  obtained  directly  from  the  SEASAT  and  GEOSAT 
Geophysical  Data  Records  (GDR).  The  original  FNOC  pressure  fields  arc  63  x  63  gridded  polar  projec¬ 
tions  of  each  hemisphere  and  arc  generated  at  three  hour  intervals  for  the  northern  hemisphere  and  at 
six  hour  intervals  for  the  southern  (Parke,  1980).  Values  of  sea  level  pressure  were  interpolated  in 
space  and  time  from  the  grids  to  correspond  to  the  altimeter’s  ground  track.  The  NMC  data  used  for 
comparison  were  obtained  from  the  National  Center  for  Atmospheric  Research  in  Boulder,  Colorado. 
These  consist  of  sea  level  atmospheric  pressure,  given  twice  daily  from  September  through  October  of 
1978  and  November  through  December  of  1986,  for  a  2.5°  x  2.5°  grid  over  the  earth’s  surface. 

To  calculate  the  differences  between  these  two  data  sets,  we  partitioned  the  earth  into  a  2.5°  x  2.5° 
grid  system  like  that  of  the  NMC  data  set.  For  each  three  day  repeat  orbit,  every  pressure  value  from 
the  GDR  was  assigned  to  a  grid  cell  as  determined  by  its  location  on  the  satellite’s  ground  track.  All 
pressure  values  in  each  cell  were  then  averaged.  NMC  values,  representative  of  the  time  the  satellite 
passed  over  a  grid  cell,  were  also  averaged  over  the  three  day  time  period.  Six  74  x  145  grids  for  each 
of  the  two  data  sets  were  generated  to  represent  the  September  15*  -  17*  orbit,  the  September 
18*  -  20*  orbit,  and  the  October  5*  -  7*  orbit  of  SEASAT  and  the  November  16*  -  18*  orbit,  the 


50 


November  19“*  -  21'*  orbit,  and  the  December  1‘*  -  Yd  orbit  of  GEOSAT. 


RESULTS 


The  FNOC  and  NMC  global  pressure  grids  for  the  September  15rt  -  17'*  orbit  (Orbit  Asea)  of 
SEASAT  were  differenced  and  all  areas  with  absolute  offsets  greater  than  10,  20,  and  30  mbar  are 
shown  in  Figure  V.l.  Figure  V.2  displays  the  offsets  between  these  two  data  sets  for  the  November 
16"*  —18**  orbit  (Orbit  AaE0)  of  GEOSAT  (Note  that  for  the  GEOSAT  data  the  increments  arc  in  units 
of  5  mbar  instead  of  the  10  mbar  increments  used  for  the  SEASAT  data).  Large  discrepancies  exist 
over  much  of  the  southern  ocean  with  differences  of  up  to  34  and  16  mbar  concentrated  in  the  Southern 
Pacific  for  the  SEASAT  and  GEOSAT  missions  respectively. 

We  wondered  whether  the  discrepancies  were  time  independent  or  not.  To  study  this  question,  we 
used  the  NMC  and  FNOC  pressure  grids  for  the  second  orbit,  (September  18'*-20'*,  Orbit  BSEa  ,  for 
SEASAT  and  November  19**— 21'* ,  Orbit  B GE0 ,  for  GEOSAT)  and  differenced  them  with  the 
corresponding  grids  from  Orbit  Asea  and  Orbit  AGEO  respectively.  Again,  the  resulting  NMC  and 
FNOC  grids  were  compared.  The  results,  shown  in  Figure  V.3  and  Figure  V.4,  thus  represent 
differences  between  the  FNOC  and  NMC  representations  of  the  variability  in  pressure  over  only  3  days. 
There  are  many  instances  where  the  differences  are  greater  than  10  mbar.  The  largest  offset  in  the 
SEASAT  data,  33  mbar,  and  the  GEOSAT  data,  21  mbar,  are  again  located  in  the  Southern  Pacific. 

An  even  larger  discrepancy  was  found  by  subtracting  the  Orbit  A  SEA  pressure  from  the  pressure  fields 
for  the  October  5'* -8'*  orbit  (Orbit  CSkA).  The  resulting  differences  between  the  NMC  and  FNOC  in 
that  case  reached  41  mbar  in  the  Southern  Pacific  (Figure  V.5).  This  result  reflects  differences  in  the 
predicted  variability  over  roughly  three  weeks. 

A  differencing  of  the  pressure  fields  for  the  December  l'*-3r<i  orbit  of  GEOSAT  (Orbit  CGE0 )  and 
Orbit  Ageo  is  shown  in  Figure  V.6.  The  differences  in  this  case  reached  almost  20  mbar  in  the  South¬ 
ern  Pacific. 


CONCLUSIONS 


These  results  demonstrate  large  differences  between  the  FNOC  and  NMC  pressure  fields  over  the 
southern  ocean  during  the  Fall  of  1978  and  the  Winter  of  1986,  with  the  biggest  discrepancies  in  the 
Southern  Pacific.  The  differences,  which  arc  undoubtably  related  to  the  poor  spatial  coverage  of 
radiosonde  stations  in  the  southern  hemisphere  (David  Salstcin,  personal  communication),  are  strongly 
lime  dependent.  For  example,  the  FNOC  and  NMC  data  sets  were  found  to  disagree  by  more  than  30 
mbar  in  their  representation  of  pressure  variability  over  only  three  days  during  the  SEASAT  experiment. 
Errors  of  this  magnitude  in  the  sea  level  pressure  fields  can  cause  errors  of  more  than  30  cm  when 
removing  pressure-driven  changes  in  sea  level  from  satellite  altimeter  data.  This  can  severely  degrade 
altimeter  estimates  of  wind-driven  circulation  in  die  southern  ocean.  Analysis  of  the  more  recent 
GEOSAT  data  indicate  better  agreement  between  the  two  data  sets  but  the  offsets  are  still  large  enough 
to  warrant  concern.  These  results  suggest  the  need  for  a  careful  assessment  of  the  sea  level  pressure 
fields  provided  by  the  various  meteorological  centers  around  the  world,  with  emphasis  on  how  those 
fields  can  be  improved  over  the  southern  ocean. 


CHAPTER  V  REFERENCES 


Parke,  M.,  G.H.  Bom,  J.F.  Scott,  1980;  SEASAT:  Altimeter  Geophysical  Algorithm  Specifications,  JPL 
Tech.  Pub.  # 622-226 . 


51 


CHAPTER  V  FIGURE  CAPTIONS 


Figure  V.l:  Differences  in  NMC  and  FNOC  sea  level  pressure  values  for  a  time  period  corresponding  to 
the  first  three  day  repeat  orbit  of  SEAS  AT  (September  15-17,  1978).  The  largest  differences  are  con¬ 
centrated  in  the  Southern  Pacific  ocean  and  reach  -34  mbar. 

Figure  V.2:  Differences  in  NMC  and  FNOC  sea  level  pressure  values  for  a  time  period  corresponding  to 
the  first  three  day  repeat  orbit  of  GEOS  AT  (November  16-18,  1986).  The  largest  differences  are  con¬ 
centrated  in  the  Southern  Pacific  ocean  and  reach  16  mbar. 

Figure  V.3:  Time  dependent  differences  in  the  NMC  and  FNOC  sea  level  pressure  values.  Data  points 
were  obtained  by  removing  the  pressure  values  of  the  first  three  day  repeat  orbit  of  SEASAT  (Sep¬ 
tember  15-17,  1978)  from  the  second  three  day  orbit  (September  18-20).  The  largest  value  of  the 
data  plotted  here  is  33  mbar. 

Figure  V.4:  Time  dependent  differences  in  the  NMC  and  FNOC  sea  level  pressure  values.  Data  points 
were  obtained  by  removing  the  pressure  values  of  the  first  three  day  repeat  orbit  of  GEOSAT 
(November  16-18,  1978)  from  the  second  three  day  orbit  (November  19-21).  The  largest  value  of  the 
data  plotted  here  is  24  mbar. 

Figure  V.5:  Same  as  Figure  V.3  but  with  the  first  orbit  having  been  removed  from  a  later  three  day  orbit 
(October  5-7,  1978)  in  the  mission.  The  largest  difference  is  -41  mbar. 

Figure  V.6:  Same  as  Figure  V.4  but  with  the  first  orbit  having  been  removed  from  a  later  three  day  orbit 
(December  1-3,  1978)  in  the  GEOSAT  mission.  The  largest  difference  is  -41  mbar.  *  Note  that  for 
Figure  V.2,  Figure  V.4  and  Figure  V.6  the  divisions  arc  in  increments  of  5  instead  of  10 
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Figure  V.2 
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CHAPTER  VI 


SPECTROSCOPIC  ANALYSIS  OF  GLOBAL  TIDE  GAUGE  SEA  LEVEL  DATA 


SUMMARY 


Yearly  and  monthly  tide  gauge  sea  level  data  from  around  the  globe  are  fitted  to  numerically  gen¬ 
erated  tidal  data  to  search  for  the  18.6  year  lunar  nodal  tide  and  14  month  pole  tide.  Both  tides  are 
clearly  evident  in  the  results,  with  amplitudes  and  phases  that  are  consistent  with  a  global  equilibrium 
response.  Global  atmospheric  pressure  data  arc  fitted  to  global,  monthly  sea  level  data  to  study  the 
response  of  the  ocean  to  pressure  fluctuations.  The  global  response  of  sea  level  to  pressure  is  found  to 
be  inverted  barometer  at  periods  greater  than  two  months. 

A  large  coherence  at  437  days  between  pressure  and  sea  level  data  is  found  for  the  North  Sea,  Baltic 
Sea  and  the  Gulf  of  Bothnia.  The  results  arc  not  entirely  conclusive,  but  they  tend  to  support 
O’Connor’s  suggestion  [1986]  that  the  apparent  enhanced  pole  tide  in  these  basins  may  be  due  to 
meteorological  forcing  rather  than  to  a  basin-scale  resonance. 

Finally,  global  averages  of  tide  gauge  data,  after  correcting  for  the  effects  of  post  glacial  rebound  on 
individual  station  records,  reveal  an  increase  in  sea  level  over  the  last  80  years  of  between  1.1  mm/yr 
and  1.9  mm/yr.  As  part  of  the  process  of  removing  the  effects  of  post-glacial  rebound,  we  fit  those 
effects  to  the  global  tide  gauge  data  and  obtain  very  good  agreement  with  results  predicted  from  post¬ 
glacial  rebound  models.  This  tends  to  support  the  post-glacial  model  results  and  it  suggests  that  the 
global  tide  gauge  data  are,  indeed,  capable  of  resolving  changes  in  sea  level  at  the  mm/yr  level. 


INTRODUCTION 


Long  period  ocean  tides  affect  estimates  of  certain  geophysical  parameters,  in  some  cases  through 
oceanic  contributions  to  the  Earth’s  inertia  tensor  and,  in  others,  through  crustal  deformation  caused  by 
the  tidal  loading.  Two  examples  of  such  estimates  arc  the  use  of  satellite  solutions  for  the  earth’s  J2 
gravity  coefficient  to  constrain  mantle  anelasticity  at  the  lunar  tidal  period  of  18.6  years,  and  the  use  of 
the  observed  period  and  damping  of  the  Chandler  wobble  to  estimate  anelasticity  at  the  14  month  wob¬ 
ble  period. 

The  18.6  year  solid-Earth  and  ocean  tides,  related  to  the  precession  of  the  lunar  nodes,  cause  an  18.6 
year  variation  in  J2.  Lambcck  and  Nakiboglu  [1983]  assumed  the  ocean  tide  was  equilibrium  (for  an 
’equilibrium’  tide,  the  sea  surface  is  a  surface  of  constant  potential)  and  used  Rubincam’s  [1984] 
observed  18.6  year  variability  in  J2  to  constrain  the  value  of  mantle  anelasticity  at  this  period.  The 
result  was  surprisingly  large.  To  help  assess  the  importance  of  the  equilibrium  ocean  tide  assumption 
on  the  conclusion,  we  note  that  if  the  earth  were  assumed  to  be  elastic,  Rubincam’s  result  would  imply 
an  18.6  year  ocean  tide  that  was  twice  equilibrium.  The  J2  result  should  continue  to  improve  as  more 
satellite  data  are  acquired. 

The  Chandler  wobble  constraint  on  anelasticity  [see,  for  example,  Smith  and  Dahlen,  1981]  depends 
critically  on  the  response  of  the  ocean  to  the  incremental  centrifugal  force  caused  by  the  wobble.  That 
response,  called  the  pole  tide,  is  known  to  affect  the  period  and,  if  it  were  non-equilibrium,  could  con¬ 
tribute  to  the  damping.  In  fact,  a  global  departure  from  equilibrium  of  only  10  percent  could  have 
observable  consequences.  Theoretical  models  for  the  pole  tide  in  the  deep  ocean  suggest  the  pole  tide 
should  be  equilibrium  [O’Connor  and  Starr,  1983;  Carton  and  Wahr,  1986;  Dickman,  1988a].  However, 
the  pole  tide  in  the  North  and  Baltic  Seas  and  in  the  Gulf  of  Bothnia  is  apparently  several  times  the 
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equilibrium  value  [Miller  and  Wunsch,  1973].  Wunsch  [1986]  postulated  that  this  enhanced  tide  may 
be  due  to  a  local  oceanic  resonance  in  the  region.  Alternatively,  O’Connor  [1986]  suggested  that  the 
apparent  14  month  spectral  peak  could  be  caused  by  meteorological  forcing  and  is  not  related  to  the 
pole  tide. 

The  response  of  the  ocean  to  variations  in  atmospheric  pressure  is  another  oceanographic  disturbance 
with  implications  for  solid  earth  geophysics.  An  inverted  barometer  response  (ic.  a  10  mm  depression 
of  sea  level  for  every  1  mbar  increase  in  pressure)  is  suggested  by  analytical  models  [see,  eg.,  Wunsch, 
1972;  Munk  and  MacDonald,  1975;  Dickman,  1988b].  The  inverted  barometer  assumption  has  been 
invoked  when  studying  the  effects  of  atmospheric  pressure  on  the  Earth’s  rotation  [see,  eg.,  Munk  and 
Hassan,  1961;  Wilson  and  Haubrich,  1976;  Merriam,  1982;  Wahr,  1983]  and  when  estimating  the  cru¬ 
stal  deformation  caused  by  pressure  fluctuations  [Rabbel  and  Zschau,  1985;  Van  Dam  and  Wahr,  1987], 
The  results  in  both  cases  are  sensitive  to  the  accuracy  of  this  assumption.  It  is  also  important  to  under¬ 
stand  the  pressure  response  in  order  to  interpret  observed  sea  level  variability  in  terms  of  wind-driven 
ocean  currents.  For  example,  oceanographers  who  process  altimeter  data  to  learn  about  the  wind-driven 
response,  usually  remove  an  inverted  barometer  pressure  response  during  data  pre-processing. 

Climate  models  that  are  used  to  study  the  effects  of  atmospheric  greenhouse  gasses  predict  an 
increase  in  the  global  temperature  over  the  next  century  of  from  1  to  4  degrees  centigrade  [Hansen,  et 
al.,  1981].  An  increase  of  this  magnitude  could  have  numerous  catastrophic  effects,  not  the  least  of 
which  would  be  a  global  rise  in  sea  level  due  to  a  combination  of  melting  polar  ice  caps  and  the  ther 
mal  expansion  of  sea  water.  One  of  the  important  goals  of  global  change  studies  is  to  improve  our 
understanding  of  this  variability  in  water  storage. 

There  have  been  previous  attempts  to  use  tide  gauge  data  to  constrain  the  18.6  year  and  14  month 
tides  [see,  eg.,  Munk  and  Cartwright,  1966;  Cartwright  and  Taylcr,  1971;  Miller  and  Wunsch,  1973;  aid 
Hosoyama,  et  al.,  1976]  the  response  of  the  ocean  to  atmospheric  pressure  [Chelton  and  Enfield  1986] 
and  the  global  rise  in  sea  level  [Emery,  1980;  Gomitz  and  Lebcdcff,  1982;  Bameu,  1983;  Peltier, 
1986].  These  studies  have  primarily  involved  the  analysis  of  data  from  individual  tide  gauges  or,  at 
most,  from  a  small  subset  of  all  available  stations. 

In  the  present  study,  we  combine  tide  gauge  data  from  several  hundred  stations  scattered  around  the 
globe,  to  test  the  hypotheses  that  the  global  18.6  year  and  14  month  tides  arc  equilibrium,  and  that  the 
response  to  pressure  is  inverted  barometer.  Our  tide  gauge  data  set  consists  of  monthly  sea  level  values 
obtained  from  the  Permanent  Service  for  Mean  Sea-Level  (PSMSL),  at  Bidston,  England  [see  Pugh  and 
Faull,  1983a].  The  distribution  of  the  721  stations  over  the  globe  is  shown  in  [Pugh  and  Faull,  1983b]. 
In  all  cases,  our  observational  results  arc  consistent  with  these  hypotheses.  (Although  the  time  resolu¬ 
tion  of  the  data  limits  the  investigation  of  the  pressure-driven  response  to  forcing  periods  of  greater  than 
two  months). 

We  investigate  the  possible  causes  of  the  apparently  enhanced  pole  tide  in  the  North  Sea  region.  Our 
comparisons  between  tide  gauge  and  pressure  data  tend  to  support  the  hypothesis  that  the  spectral  peak 
io  due  to  narrow-band  meteorological  lorcing  in  the  region,  rather  than  to  a  basin-scalc  resonance.  Our 
conclusions,  though,  should  be  regarded  as  tentative,  because  the  dominant  meteorological  forcing 
would  be  due  to  winds,  but  we  did  not  include  wind  data  in  our  analysis. 

We  also  use  the  PSMSL  data  to  constrain  the  global  rise  in  sea  level.  Although  our  estimates  vary 
somewhat,  depending  on  how  we  restrict  the  data  and  correct  for  post  glacial  rebound  [sec  eg.,  Peltier, 
1986]  we  infer  a  global  sea  level  rise  of  between  1.1  mm/yr  and  1.9  mm/yr,  with  a  preferred  value  of 
1.75  mm/yr  for  stations  in  the  regions  covered  by  Peltier’s  model.  In  the  process  of  correcting  for  the 
effects  of  post  glacial  rebound,  we  least  squares  fit  Peltier’s  [1986]  predicted  vertical  uplift  to  the  data, 
and  find  very  good  agreement  (ie.  a  fit  with  a  value  close  to  1). 


METHODS  OF  ANALYSES 


In  this  section,  we  describe  two  methods  wc  use  to  help  identify  small  signals  in  the  global  tide 
gauge  data.  Wc  use  the  least  squares  fit  technique  to  study  the  ocean’s  response  to  atmospheric  pres¬ 
sure.  And,  we  stack  the  data  to  investigate  the  18.6  year  tide  and  the  14  month  tide.  We  use  both  tech¬ 
niques  to  study  the  global  rise  in  sea  level  and  the  effects  of  post  glacial  rebound.  We  use  the  monthly 
PSMSL  data  to  study  the  pole  tide  and  the  response  to  atmospheric  pressure.  We  use  yearly  averages 
of  the  PSMSL  data  to  study  the  18.6  year  tide  and  the  global  rise  in  sea  level. 


Least  squares  fit 

Let  s*(r)  and  pk{i)  represent  the  time  dependent  fluctuations  in  sea  level  and  in  pressure  at  tide 
gauge  k .  Separate  the  observed  sea  level  fluctuations  into  a  component  caused  by  the  pressure,  and  a 
remainder,  e*(r),  due  to  a  combination  of  wind-driven  fluctuations,  tides,  and  observational  noise.  Sup¬ 
pose  that  the  pressure-induced  variability  at  an  individual  tide  gauge  depends  on  the  past,  present,  and 
future  values  of  pressure  as  measured  at  that  gauge  alone.  (Except  in  cases  where  the  arrival  of 
pressure-driven  waves  precedes  the  arrival  of  the  change  in  pressure,  causality  precludes  dependence  on 
the  future,  but  by  solving  for  that  dependence  we  can  partially  assess  the  accuracy  of  our  results.)  Then, 
the  most  general  linear  relationship  between  sk  and  pk  is: 


•s*(0  =  |A{t)p*(/-t)  dx  +  £*(/) 


(1) 


Here,  A  (x)  can  be  interpreted  physically  as  a  Green’s  function  describing  the  response  of  the  ocean  to 
an  impulsive  pressure  disturbance  (the  true  Green’s  function  would  actually  be  more  complicated  than 
A(x),  since  equation  (1)  assumes  that  sea  level  at  any  point  is  unaffected  by  pressure  disturbances  at 
other  points).  Since  the  ocean  has  a  finite  memory,  A  (t)  should  approach  zero  for  large  t. 

Suppose  that  both  sk(t)  and  pk(i)  are  discretized  to  a  time  series  of  monthly  values,  i  =f, ,  and  sup¬ 
pose  A  (x)  is  negligible  for  lit  >  Lf, ,  where  L  is  an  integer.  If  At  is  the  sampling  interval  (equal  to  one 
month  in  our  case),  then  (1)  reduces  to  the  discrete  form: 

L 

h ( tj )  =  E  A  (/,  )pk ( tj -u )A t  +ek(tj)  (2) 


For  an  inverted  barometer  ocean,  A(0)  would  be  -10.1  mm/mbar  [see  Wunsch,  1972],  and  all  other 
A  (tt )  would  be  zero. 

Suppose  A(r.)  is  independent  of  the  tide  gauge  location.  This  is  equivalent  to  assuming  that  the 
ocean  responds  to  pressure  in  the  same  way  at  every  location.  Then,  we  can  estimate  A  (t,)  for  each 
by  least  squares  fitting  to  all  PSMSL  tide  gauge  data  simultaneously.  We  can  least  squares  fit  each 
A  (tj )  individually,  or  all  the  A  (t, )  simultaneously.  As  an  example,  if  each  data  point  is  given  equal 
weight,  regardless  of  where  or  when  it  was  taken,  and  if  Mk  is  the  total  number  of  months  in  the  time 
series  for  tide  gauge  k ,  then  the  least  squares  fit  for  an  individual  A  (t; )  gives: 

E  E  Sk(tj)Pk(tj-ti) 

AM  =  — -  (3) 

E  E 

*=g«uge  /«! 
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As  another  application,  the  PSMSL  data  arc  simultaneously  fit  to  theoretically  predicted  post  glacial 
rebound  data  and  to  a  global  linear  trend,  to  estimate  the  global  rise  in  sea  level  and  to  help  assess  the 
rebound  models.  The  rebound  data  are  constructed  from  the  linear  trends  predicted  for  each  station  by 
the  models,  and  are  identical  to  the  PSMSL  data  in  coverage  and  gaps.  Whereas  the  effects  of  ice 
melting  and  thermal  expansion  on  sea  level  is  assumed  here  to  be  independent  of  position,  the  linear 
trends  due  to  post  glacial  rebound  are  not.  To  the  extent  the  rebound  varies  with  position,  it  may  be 
separated  from  the  general  increase  in  sea  level  by  the  simultaneous  fit  The  fit  to  the  rebound  model 
gives  information  on  which  stations  should  be  included  in  the  simultaneous  fit,  as  described  below.  It 
also  allows  us  to  partially  assess  the  accuracy  of  the  rebound  models,  and  of  the  global  sea  level  rise 
inferred  from  the  tide  gauge  data. 


Slacking 

The  procedure  of  stacking  multi-station  data  to  enhance  small  signals  has  been  used  by  seismologists 
to  improve  their  estimates  of  the  Earth’s  free  oscillation  cigcnfrequencics  [Gilbert  and  Dzicwonski, 
1975],  and  to  search  for  short  period  oceanic  normal  modes  in  the  Pacific  [Luther,  1982].  Stacking  is 
particularly  useful  in  cases  where  the  spatial  dependence  of  the  signal  is  known  beforehand. 

Suppose  we  want  to  test  the  hypothesis  that  the  18.6  year  tide  is  equilibrium.  The  18.6  year  equili¬ 
brium  amplitude  at  co-latitude  0,  eastward  longitude  X,  and  time  /,  has  the  form  (using  the  corrected 
tables  of  tidal  harmonics  from  Cartwright  and  Edden  [1973]): 

H(i  ,0,X)  =  27.94  5(0, X)  cos  (tot  +  <p>  mm  (4) 

where  to  is  the  frequency  of  the  18.6-year  tide  and  d>  is  the  phase  at  r=0.  The  function  5(0  A) 
represents  the  spatial  dependence  of  the  equilibrium  tide,  and  would  equal  (l-Ht-/i)Y°  were  it  not  for 
mass  conservation  of  the  oceans,  sea  floor  loading,  and  gravitational  self  attraction  (X  and  k  arc  tidal 
Love  numbers,  and  the  spherical  harmonic  Y®  is  normalized  so  that  the  integral  over  the  unit  sphere  of 
IY2°I2  is  one).  These  latter  effects  cause  the  introduction  of  other  spherical  harmonics  into  5(0A)  and 
increase  the  Y®  component  by  about  10  percent  The  function  5(0 A)  can  be  found  by  iterative  solution 
of  equation  (102)  in  Dahlen,  [1976], 

Once  we  have  found  5(0A),  we  least  squares  fit  it  to  all  the  PSMSL  tide  gauge  data  at  each  time  t 
(see  equation  6,  below).  The  average  and  the  secular  trend  arc  removed  from  each  station  prior  to  the 
fit  The  resulting  time  series  of  fit  values  is  referred  to  here  as  a  stack.  The  stack  is  spectrally  analyzed 
to  search  for  a  peak  at  the  18.6  year  frequency,  to.  If  a  peak  is  found,  its  amplitude  and  phase  can  be 
compared  with  the  results  that  would  be  expected  for  an  equilibrium  tide.  To  estimate  these  equili¬ 
brium  results,  an  synthetic  data  set  is  constructed  by  replacing  the  PSMSL  sea  level  value  at  every 
gauge  and  at  every  time  by  the  estimated  equilibrium  value,  h(t,9,X)  in  (4),  at  that  same  position  and 
time.  The  synthetic  data,  like  the  PSMSL  data,  are  stacked  against  the  equilibrium  spatial  dependence 
5(0,X)  and  the  amplitude  and  phase  of  the  spectral  peak  at  18.6  years  arc  estimated. 

Good  agreement  between  these  stacks  of  PSMSL  and  synthetic  data  is  not  enough  to  guarantee  that 
the  tide  is  equilibrium.  It  is  desirable  to  also  stack  against  other  functions  of  0  and  X.  For  example, 
the  equilibrium  5(9A)  for  the  18.6  year  tide  is  nearly  proportional  to  Y®(0,X).  Wc  stack  against  other 
spherical  harmonics  as  well,  up  to  degree  6.  These  additional  stacks  serve  two  purposes.  First,  by 
stacking  the  synthetic  data  we  can  assess  the  effectiveness  of  the  stacking  procedure.  For  instance,  is 
the  station  coverage  complete  enough  that  different  spherical  harmonics  arc  reasonably  orthogonal?  Wc 
find  that  it  is.  Second,  by  stacking  the  PSMSL  data  and  looking  at  the  spectra,  we  can  further  test  the 
equilibrium  assumption.  For  example,  if  the  data  are  stacked  against  a  spheri'-M  harmonic  other  than 
Y®,  and  if  the  18.6  year  tide  is  equilibrium,  then  there  should  not  be  a  notable  ^^ak  at  18.6  years  in  the 
spectrum  of  the  stack.  Although  we  do  not  show  results  here  for  these  other  stacks,  wc  find  no  evi¬ 
dence  of  18.6  year  peaks  that  stand  significantly  above  the  noise  continuum  for  any  stacks  against  pure 
spherical  harmonics  other  than  Y® .  Wc  compare  these  PSMSL  results  with  Y"  stacks  of  the  synthetic 
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data,  and  find  that  the  18.6  year  peaks  in  the  synthetic  slacks  are  also  below  the  noise  continuum  of  the 
PSMSL  spectra.  These  results,  then,  arc  consistent  with  the  equilibrium  hypothesis.  That  hypothesis  is 
addressed  further,  below,  when  we  discuss  the  stacks  against  the  equilibrium  C(9,X). 

A  similar  procedure  is  used  to  investigate  the  14  month  pole  tide.  The  equilibrium  pole  tide  is  gen¬ 
erated  from  polar  motion  data  (we  use  polar  motion  values  from  the  International  Latitude  Observatory 
of  Mizusawa;  see  Yumi  and  Yokoyama,  1980)  according  to: 

//(/  ,8,X)  =  Rc  [[«,(!)  +  im2(o]c(0,X.)]  (5) 

Here,  mt(i)  and  m2(t )  are  the  x  and  y  coordinates  of  the  pole  position  (with  magnitude  up  to  one  third 
of  an  arc  second),  Q  and  a  arc  the  Earth’s  mean  rotation  rate  and  radius,  and  Re  denotes  the  real  part. 
The  spatial  dependence,  £(0,X),  in  (5)  is  complex,  and  so  when  it  is  least  squares  fit  to  the  PSMSL  data 
it  gives  a  complex  stack.  In  the  absence  of  sea  floor  loading,  gravitational  self  attraction,  and  mass 
conservation,  C(0,X)  would  equal  (1  +k-h)Y2(QX).  As  in  the  case  of  the  18.6  year  tide,  is 

modified  by  about  10%  due  to  these  additional  effects,  and  it  can  be  estimated  as  described  by  Dahlcn 
(1976).  Stacks  against  the  equilibrium  £(9,X)  arc  described  below.  Stacks  of  the  PSMSL  data  against 
spherical  harmonics  other  than  K2  do  not  exhibit  significant  peaks  at  the  14  month  pole  tide  period,  and 
are  consistent  with  Y"  slacks  of  the  synthetic  data. 

For  any  of  these  examples,  the  least  squares  fit  to  form  a  stack  can  either  be  weighted  or  unweighted. 
One  weighting  procedure  we  find  to  be  particularly  useful  involves  stacking  on  a  grid.  The  earth's  sur 
face  is  divided  into  grid  elements,  in  our  case  either  10°  latitude  by  20°  longitude  or  15°  latitude  by  30° 
longitude.  Sea  level  heights  for  each  station  (with  secular  trends  and  station  averages  removed)  are 
multiplied  by  the  equilibrium  spatial  dependence,  £(0,X),  and  the  results  are  averaged  over  individual 
grid  elements  for  a  given  time,  t.  A  global  average  is  then  constructed  by  summing  the  grid  averages. 
The  final  result  is  normalized  by  dividing  by  the  global  average  of  the  equilibrium  data.  The  resulting 
time  scries  is  the  weighted  stack. 

For  example,  let  sk  (t )  be  the  tide  gauge  reading  for  station  k  at  time  t,  with  the  secular  trend  and  sta¬ 
tion  average  removed.  Let  0t  and  Xk  be  the  co-latitude  and  eastward  longitude  of  station  k.  Let  np(i) 
be  the  total  number  of  gauge  measurements  in  grid  element  p  at  time  t ,  and  suppose  there  are  N  grid 
points.  Then,  the  complex  gridded  stack,  x(t ),  is  defined  as 

X-777  X  hiOWM 

p=\  np\*  )  *=gauges 

*(D  =  -W~. - — -  (6) 

Z—T  X  C(0*.^)C(0*A) 

p=\  ™p\*)  tsgauges 
in  p 


Here  *  denotes  complex  conjugation.  This  result,  (6),  is  equivalent  to  a  weighted  least  squares  fit  of  £ 
to  the  data  at  time  t ,  and  it  reduces  to  an  ungridded,  unweighted  stack  when  the  grid  elements  are  so 
small  that  np(t)  =  1  for  every  p. 

Individual  station  records  exhibit  considerable  long  period  variability,  most  of  it  probably  caused  by 
forcing  from  surface  winds.  This  makes  it  difficult  to  use  a  single  tide  gauge  record  to  identify  any  of 
the  long  wavelength  spectral  features  we  arc  considering  here.  Stacking  the  individual  stations  against 
the  spatial  dependence  of  the  equilibrium  solution  amplifies  any  signal  with  that  spatial  dependence,  and 
this  can  significantly  reduce  the  correlation  with  the  wind-driven  variability. 

The  use  of  a  grid  de-emphasizes  areas  with  great  concentrations  of  stations,  such  as  northern  Europe 
and  Japan,  by  weighting  them  equally  with  areas  of  fewer  stations.  Gridding  also  tends  to  further  reduce 
the  contamination  from  long  wavelength,  wind-driven  features.  Roughly  speaking,  the  contamination  of 
the  final  stack  is  reduced  when  the  wind-driven  signal  has  wavelengths  much  longer  than  the  grid  size. 
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but  is  increased  at  shorter  wavelengths. 

For  example,  imagine  a  situation  where  there  are  10  stations:  9  within  grid  element  A  and  1  within  a 
far  distant  grid  element  B.  If  the  wind-driven  signal  has  a  wavelength  much  longer  than  the  grid  size, 
then  within  grid  element  A  that  signal  is  highly  correlated  with  the  long  wavelength  feature  we  are 
interested  in:  both  are  nearly  constant  across  the  element.  In  that  case,  an  ungridded  stack  would  not 
effectively  remove  the  wind-driven  signal,  since  the  stack  would  be  dominated  by  the  high  correlation 
at  the  9  stations  within  A.  A  gridded  stack  would  be  more  effective,  since  then  the  9  elements  in  A 
would  be  weighted  equally  with  the  1  element  in  B,  and  thus  the  stack  could  more  effectively  exploit 
any  differences  between  the  signals  at  A  and  B. 

On  the  other  hand,  since  isolated  stations  contribute  more  heavily  to  a  gridded  stack,  that  stack  is 
more  sensitive  to  noise  in  those  individual  station  records.  The  sizes  of  the  grids  used  here  are  chosen 
to  be  large  enough  to  minimize  the  effects  of  wind-driven  fluctuations  and  of  gaps  in  the  grid  averages, 
but  smaller  than  the  wavelength  of  the  dominant  spherical  harmonic  in  the  stacking  function.  The 
yearly  data  are  stacked  on  a  coarser  grid  than  are  the  monthly  data.  In  each  case,  the  grid  size  that  is 
selected  is  the  one  for  which  the  background  noise  for  the  stack  is  minimized. 

Subtracting  the  average  tide  gauge  height  for  each  station  prior  to  a  stack  is  necessary  to  establish  a 
uniform  benchmark  for  all  stations.  But,  if  a  station  docs  not  report  data  over  the  entire  time  span  of 
the  stack,  the  computed  station  average  is  not  the  true  average,  and  this  can  introduce  systematic  errors 
into  the  spectrum.  (Of  the  721  stations  in  the  PSMSL  data  set,  only  7  cover  the  entire  interval  between 
1900  and  1979.)  The  frequency  content  of  these  errors  depends  on  the  time  spans  of  the  individual  sta¬ 
tion  data.  The  errors  are  likely  to  be  important  at  periods  close  to  and  longer  than  the  average  of  the 
station  lengths.  Although  we  have  found  no  consistent  way  of  eliminating  these  errors,  we  have 
attempted  to  minimize  their  impact  by  only  including  stations  with  data  spans  that  arc  longer  than  the 
period  of  the  signal  wc  arc  looking  for  (at  least  when  we  can  define  that  period).  The  pole  tide  results 
are  not  particularly  sensitive  to  this  problem  since  the  average  station  length  for  all  PSMSL  stations  is 
29  years,  which  is  many  times  the  period  of  the  tide. 

Fitting  and  removing  secular  trends  from  stations  prior  to  stacking  is  found  to  remove  another  source 
of  error,  as  power  associated  with  these  trends  can  leak  into  the  spectral  band  of  interest.  In  addition, 
before  any  spectral  analysis  is  performed  on  the  final  stacks,  a  secular  trend  and  an  average  are 
removed  from  the  final  time  scries.  (Even  if  trends  and  averages  are  removed  from  the  individual  sta¬ 
tion  data,  there  can  still  be  a  trend  and  a  non-zero  average  in  the  stacked  data,  if  the  stack  is  weighted, 
because  an  individual  station  docs  not,  in  general,  report  data  over  the  entire  time  span  of  the  stack.) 

The  spectral  methods  employed  here  to  analyze  a  stack  include  finding  conventional  amplitude  and 
coherence  spectra,  and  using  the  more  sophisticated  multi-taper  method  developed  by  Thomson  [1982] 
to  minimize  spectral  leakage  and  adapted  for  geophysical  applications  by  Lindbergh  ct  al.  [1987],  For 
the  multi-taper  technique,  the  data  arc  multiplied  by  one  of  six  prolate  spheroidal  sequences,  or  tapers, 
to  create  six  lime  series.  A  spectrum  is  then  constructed  from  the  discrete  Fourier  transforms  of  each  of 
these  scries  following  section  3,  part  1,  of  Lindbergh,  cl  al.  [1987],  A  statistical  F-tcst  provides  a 
confidence  estimate  that  an  apparent  periodicity  in  the  spectrum  is  truly  sinusoidal. 

The  conventional  amplitude  spectra  shown  here  all  employ  a  cosine  taper.  Following  convention,  wc 
construct  complex  Fourier  transforms  as  twice  the  real  transform,  and  the  amplitude  spectrum  is  the 
absolute  value  of  the  complex  Fourier  transform.  Al!  of  the  amplitude  spectra  shown  in  the  figures  arc 
ovcrrcsolvcd.  By  this,  wc  mean  that  the  absolute  amplitude  is  computed  at  more  frequent  intervals  than 
the  elementary  bandwidths,  in  order  to  obtain  a  clearer  idea  of  where  the  actual  peaks  lie.  To  find  the 
coherence  at  frequency  w  between  two  lime  scries,  the  time  scries  arc  multiplied  by  a  cosine  taper,  and 
the  average  coherence  across  seven  elementary  bandwidths  centered  on  co  is  calculated  according  to 
(equation  (5)  of  Wilson  and  Haubrich  [1976]).  Wiui  this  choice  for  the  coherence  spectrum,  wc  can  be 
95%  confident  that  two  scries  arc  correlated  al  to,  if  the  coherence  at  0)  is  larger  than  0.43  .  Other 
choices  of  tapers  and  bandwidths  arc  possible  [sec  Chao,  1988]. 
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RESULTS 


The  response  to  pressure. 


Figure  VI. 1  shows  a  plot  of  the  least  squares  fit  parameter,  A  (z)  (equation  (3)  with  t,=z  ),  between 
sea  level  and  atmospheric  pressure,  for  lag  times  of  up  to  ±  240  months.  The  pressure  data  are  obtained 
from  an  objective  analysis  of  global  station  pressure  data  for  January,  1900  -  April,  1973  [Wahr,  1983], 
and  so  we  restrict  the  PSMSL  data  to  this  time  period.  Seasonal  effects  are  removed  from  both  sea 
level  and  pressure  data  sets  prior  to  the  fit. 

All  PSMSL  stations  during  this  time  period  arc  included  in  the  fit  shown  in  Figure  VI.  1,  except  for 
those  in  the  North  and  Baltic  Seas  and  in  the  Gulf  of  Bothnia.  As  shown  below,  the  Baltic  Sea  and  Gulf 
of  Bothnia  do  not  appear  to  be  inverted  barometer  (although  the  North  Sea  does).  Furthermore,  by 
omitting  these  3  shallow  seas  from  the  data,  we  end  up  using  the  same  data  subset  in  Figure  VI.  1  that 
we  will  use  later  to  study  the  global  pole  tide. 

The  result  for  A  (0)  is  -10.1  mm/mbar  for  the  fit  shown  in  Figure  VI. 1.  The  results  for  all  other  A  (t) 
are  close  to  zero,  with  an  rms  value  of  .5  mm/mbar.  Consequently,  we  assign  ±  .5  mm/mbar  to  be  the 
uncertainty  in  our  estimate  of  A  (0).  Since  wind  data  were  not  analyzed,  there  is  no  way  for  us  to 
separate  the  contribution  of  winds  to  the  fit  to  pressure.  There  should  be  a  temporal  correlation 
between  winds  and  pressure  at  each  each  tide  eauge.  Our  hope  is  that  there  are  enough  stations  over  a 
wide  enough  area  in  the  global  average  to  break  the  spatial  correlation  between  winds  and  pressure. 
Without  wind  data  we  cannot  confirm  that  this  has  happened,  but  the  results  shown  in  Figure  VI. 1  are 
at  least  consistent  with  an  inverted  barometer  response  of  the  oceans  to  pressure  at  periods  of  two 
months  and  longer. 

The  response  for  zero  lag  for  other  sets  of  stations  is  as  follows: 

center,  delim  0.  tab(');  l  1  c.  Global  oceans  including  all  shallow  basins:! A  (0)=  -12.1  ±  .9 
mm/mbar,!574  stations  Northern  hemisphere  excluding  North  Sea  region:! A  (0)=  -10.0  ±  .5 
mm/mbar,!407  stations  Eastern  hemisphere  excluding  North  Sea  region:!A  (0)=  -10.3  ±  .7  mm/mbar,!234 
stations  Southern  hemisphere: !A(0)=  -10.9  ±  1.1  mm/mbar,!80  stations  Western  hemisphere: ! A  (0)=  -9.9 
±  .6  mm/mbar,!253  stations  North  Sca:!A(0)=  -9.7  ±  .9  mm/mbar,!26  stations  Baltic  Sca:!A(0)=  -12.5  ± 
1.3  mm/mbar,!42  stations  Gulf  of  Bothnia: ! A  (0)=  -16.2  ±  1.7  mm/mbar,!  19  stations  Mediterranean 
Sea:!A(0)=  -14.4  ±  1.2  mm/mbar,!33  stations 


In  all  these  fits,  the  A  (x)  for  t*0  are  close  to  zero.  The  rms  values  for  these  other  A  (t)  are  shown  as 
the  uncertainties  for  A  (0).  Note  that  as  the  number  of  stations  included  in  the  fit  decreases,  the  rms 
value  for  A  (t)  with  t*0  increases. 

These  results  include  examples  of  shallow  basins  where  the  deviation  from  the  inverted  barometer 
response  for  A  (0)  exceeds  the  rms  value  of  the  other  A  (x)’s.  This  docs  not  necessarily  imply  a  true 
non-inverted  barometer  response,  because  the  rms  value  docs  not  adequately  represent  everything  that 
could  bias  the  fit  In  particular,  wind-driven  changes  in  sea  level  could  be  correlated  with  surface  pres¬ 
sure,  and  thus  could  affect  our  estimate  of  A  (0)  without  contributing  noticeably  to  A  (t)  for  t*0.  This 
is  apt  to  be  a  particular  problem  when  considering  small  regions  at  a  time,  such  as  the  Baltic  Sea  and 
the  Gulf  of  Bothnia,  since  then  there  is  likely  to  be  greater  spatial  coherence  between  the  wind  and 
pressure  fields. 
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The  pole  tide. 

Monthly  tide  gauge  data  are  stacked,  as  described  above,  to  study  the  pole  tide.  Before  stacking  the 
data,  we  remove  the  effects  of  atmospheric  pressure,  computed  using  the  inverted  barometer  assump¬ 
tion.  Again,  we  restrict  our  analysis  to  January,  1900  -  April,  1973.  We  find  that,  for  stations  outside 
the  North  Sea,  Baltic  Sea,  and  Gulf  of  Bothnia,  the  removal  of  pressure  affects  the  pole  tide  amplitude 
inferred  below  by  about  10  percent. 

In  this  and  the  succeeding  analysis  involving  the  pole  tide,  stations  in  the  North  and  Baltic  Seas  and 
in  the  Gulf  of  Bothnia  are  removed  from  the  data  as  there  is  an  anomalously  large  spectral  peak  for  this 
region  at  periods  close  to  the  pole  tide  period  of  434  to  437  days  [see,  for  example.  Miller  and  Wunsch, 
1973],  An  equilibrium  stack  of  26  stations  in  the  North  Sea  reveals  the  PSMSL  amplitude  at  the  437 
day  period  to  be  3.6  times  the  equilibrium  amplitude,  for  the  42  Baltic  Sea  stations  die  PSMSL  ampli¬ 
tude  is  6.2  limes  equilibrium,  and  for  the  19  stations  in  the  Gulf  of  Bothnia,  the  PSMSL  amplitude 
exceeds  the  equilibrium  amplitude  by  a  factor  of  9. 

The  beat  period  between  the  pole  tide  period  and  one  year  is  approximately  60  months.  To  minimize 
the  correlation  with  the  annual  period,  we  ch  xise  to  include  in  our  global  analysis  only  those  stations 
having  greater  than  60  months  of  data.  This  leaves  us  with  487  stations  out  of  the  original  721  in  the 
full  PSMSL  data  set. 

Figure  VI.2  (a)  shows  ungridded  amplitude  spectra  for  Fj  stacks  of  these  487  stations.  Figure  VI.2 
(b)  j,hows  the  spectra  of  the  slacks  averaged  on  a  10°  by  20°  grid.  Shown  arc  the  spectra  for  the 
PSMSL  data,  for  the  synthetic  (ie.  equilibrium)  data,  and  for  differenced  data.  The  differenced  data  are 
generated  by  subtracting  the  synthetic  time  series  from  the  PSMSL  time  scries.  Note  that  in  Figure 
VI.2  the  spectra  for  both  the  PSMSL  data  and  the  synthetic  data  show  a  double  peak  for  the  pole  tide, 
with  one  peak  at  427  days  and  the  other  at  437  days.  This  is  a  well  known  feature  of  the  polar  motion 
spectrum.  Most  of  the  signal  contributing  to  this  double  peak  is  from  data  early  in  the  time  series.  The 
two  peaks  arc  close  enough  together  that  their  resolution  requires  most  of  the  74  years  of  the  data  span. 

In  both  the  gridded  and  non-gridded  stacks  the  pole  Udc  is  clearly  evident  in  the  spectrum  for  the 
PSMSL  data  and  there  is  good  visual  agreement  with  the  synthetic  data.  The  agreement  is  particularly 
good  for  the  gridded  i  :sults.  In  both  cases,  the  differenced  data  docs  not  have  the  characteristic  double 
peak  spectrum,  indicating  that  the  PSMSL  stack  is  similar  to  the  equilibrium  stack  in  amplitude  and 
phase. 

To  evaluate  the  agreement  quantitatively,  we  compute  amplitudes  for  both  the  PSMSL  and  the  syn¬ 
thetic  data  sets.  First,  the  amplitudes  at  the  frequencies  of  the  two  apparent  pole  tide  peaks  arc  com¬ 
puted,  and  the  results  from  the  PSMSL  data  and  synthetic  data  arc  compared.  We  find,  for  the  ungra¬ 
ded  data  shown  in  Figure  VI.2  (a),  that  the  amplitude  of  the  427  day  peak  in  the  PSMSL  data  is  .99  ± 
.28  times  the  synthetic  amplitude,  and  that  of  the  437  day  peak  is  .74  ±  .27  times  the  synthetic. 

We  then  repeat  this  procedure,  but  for  the  gridded  stacks  in  Figure  V1.2  (b).  In  this  case,  the  427 
day  and  437  day  amplitude  ratios  are  1.01  ±  .18  and  .96  ±  .17  respectively.  The  uncertainties  for  both 
the  gridded  and  ungridded  results  are  obtained  from  visual  estimates  of  the  signal-to-noise  ratio  at  the 
pole  tide  frequency,  and  assuming  the  worst  case:  that  the  amplitudes  of  the  signal  and  the  noise  arc 
additive  (in  other  words,  that  the  noise  is  either  in-phase  or  180°  out-of-phase  with  the  signal).  The 
background  amplitude  used  to  obtain  the  signal -to-noisc  ratio  is  1.0  mm  for  both  the  ungridded  and 
gridded  stacks. 

To  estimate  the  phase  difference  between  the  PSMSL  and  the  synthetic  data  and  to  obtain  another 
estimate  for  the  amplitude  ratio,  we  fit  the  synthetic  data  to  the  PSMSL  data  using  a  complex  constant, 
after  a  filter  is  used  to  extract,  from  both  time  scries,  the  power  in  a  spectral  band  centered  over  the  two 
pole  tide  frequencies  (.068-. 073  c/mo).  The  power  in  the  spectral  band  of  interest  is  extracted  so  that 
the  band  at  the  annual  period  and  the  low  frequency  power  docs  not  affect  the  estimate  of  phase  (there 
is  a  residual  spectral  peak  centered  around  1  c/yr,  even  though  an  annual  term  has  been  fit  and  removed 
from  the  data).  The  phase  difference  between  the  PSMSL  and  the  synthetic  data  is  small  for  both  the 
gridded  and  ungridded  stacks.  The  PSMSL  leads  the  synthetic  data  in  the  ungridded  stacks  by  3°  ±  50° 
and  the  overall  amplitude  ratio  is  .83  ±  .27;  roughly  the  average  of  the  amplitude  ratios  taken  for  each 
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spectral  peak  separately.  For  the  gridded  data  the  PSMSL  leads  the  synthetic  data  by  0°  ±  46°  and  the 
amplitude  ratio  is  .96  ±  .18.  A  phase  lead,  8,  is  defined  here  so  that  if  the  equilibrium  tide  is  propor¬ 
tional  to  cos(aw-X),  where  X  is  eastward  longitude  and  to  is  the  (positive)  pole  tide  frequency  then  the 
observed  ocean  tide  is  proportional  to  cos((o/-X+S).  Each  degree  of  phase  lead  represents  a  time  lead 
of  approximately  1.2  days. 

To  estimate  the  uncertainty  in  the  phase,  we  note  that  the  amplitude  uncertainty  found  from  the 
signal-to  noise  ratio,  as  described  above,  mosdy  reflects  the  uncertainty  in  the  in-phase  component.  We 
assume  the  uncertainty  in  the  out-of-phase  component  is  just  as  large,  and  compute  the  resulting  bounds 
on  the  phase.  This  procedure  results  in  a  larger  uncertainty  than  that  predicted  from  the  formal  error  of 
the  fit. 

We  have  also  computed  the  multi-taper  spectrum,  as  described  above,  for  the  gridded  equilibrium 
stack.  The  results  do  not  differ  significantly  from  the  conventional  spectral  results,  probably  because 
the  cosine  taper  in  the  conventional  spectra  is  sufficient  to  exclude  the  majority  of  spectral  leakage  into 
the  pole  tide  band. 

Taken  together,  the  Yj  stacks  of  monthly  data  suggest  that  the  pole  tide  amplitude  is  .9  ±  .3  times 
the  equilibrium  amplitude,  and  that  the  actual  ocean  tide  leads  the  equilibrium  pole  tide  by  2°±48°. 


The  enhanced  pole  tide  in  the  North  Seas. 

Figure  VI. 3  (a)  shows  unlapered  amplitude  spectra  of  surface  pressure,  PSMSL  sea  level  data,  and  the 
equilibrium  pole  tide  stacked  against  Yj ,  for  the  North  Sea.  The  surface  pressure  results  are  obtained 
by  multiplying  the  pressure  data  described  above  by  the  inverted  barometer  factor  of  10  mm/mbar,  and 
then  replacing  each  tide  gauge  reading  with  the  inverted  barometer  value  at  the  time  and  location  of 
that  reading.  Figures  VI. 3  (b)  and  (c)  show  the  same  spectra  for  the  Baltic  Sea  and  Gulf  of  Bothnia 
respectively.  For  these  small  seas,  a  Yj  stack  is  not  much  different  than  a  spatial  average.  The  same 
spectra  are  plotted  for  the  487  stations  that  span  the  remainder  of  the  globe  in  Figure  VI.3  (d).  In  all 
the  figures,  amplitude  is  shown  instead  of  power  to  compare  the  three  disparate  spectra  most  clearly. 
Note  that  in  all  three  shallow  seas,  there  is  a  pronounced  sea  level  peak  in  the  pole  tide  band,  that 
stands  well  above  the  equilibrium  results. 

We  have  attempted  to  interpret  the  14  month  peak  in  the  sea  level  data  in  terms  both  of  meteorologi¬ 
cal  forcing  [O'Connor,  1986]  and  of  a  local  oceanic  resonance  [Wunsch,  1986],  We  formulate  two 
hypotheses.  One  is  that  there  is  a  resonance  that  amplifies  the  response  of  this  region  to  external  forc¬ 
ing.  Such  a  resonance  should  affect  both  the  pole  tide  and  the  response  to  meteorological  forcing.  The 
other  hypothesis  is  that  there  is  no  local  resonance,  but  that  the  spectral  peak  in  sea  level  is  due  to  a 
combination  of  an  equilibrium  pole  tide  and  a  static  response  to  meteorological  forcing  (ie.  an  inverted 
barometer  response  to  pressure  and  a  response  to  local  winds).  We  assume  the  winds  are  correlated 
with  the  pressure  in  this  frequency  band,  so  that  the  total  stacked  response  of  the  ocean  to  the  meteorol¬ 
ogy  is  proportional  to  the  stacked  pressure  data  across  this  frequency  band. 

To  assess  these  two  hypotheses,  we  first  note  that  for  the  three  shallow  basins,  there  is  a  pressure 
peak  at  the  437  day  period,  close  to  coincident  with  the  pole  tide  period.  In  fact,  the  peaks  in  the  sea 
level  spectra  appear  to  line  up  reasonably  well  with  peaks  in  the  pressure  spectra.  There  seems  to  be 
less  agreement  between  the  sea  level  and  equilibrium  pole  tide  spectra.  On  the  other  hand,  for  the  glo¬ 
bal  results  there  is  good  apparent  agreement  between  sea  level  and  the  equilibrium  pole  tide,  and  less 
similarity  between  sea  level  and  pressure. 

Coherence  results  tend  to  support  this  observation.  Coherence  spectra  between  Yj  stacks  of  sea  level 
and  pressure  are  shown  in  Figure  'T.4  (a)  for  the  North  Sea,  Figure  VIA  (b)  for  the  Baltic  Sea,  and 
Figure  VIA  (c)  for  the  Gulf  of  Bothnia.  The  reason  for  using  stacks  against  Yj  instead  of  simple  aver¬ 
ages,  is  for  consistency  with  the  global  results.  To  minimize  spectral  leakage  from  the  annual  into  the 
pole  tide  bands,  seasonal  averages  have  been  removed  from  all  data  from  which  coherence  spectra  arc 
generated.  For  all  three  shallow  basins,  there  is  a  correlation  above  the  95  per  cent  confidence  level 
between  sea  level  and  pressure  in  the  frequency  band  of  the  pole  tide  (a  coherence  estimate  greater  than 


0.43).  The  coherence  is  particularly  strong  over  the  Gulf  of  Bothnia.  Although  there  is  a  5%  chance 
that  the  large  coherence  is  spurious,  it  is  interesting  that  it  occurs  in  a  frequency  band  where  there  are 
large  amplitude  spectral  peaks  in  both  pressure  and  sea  level.  Figure  V1.4  (d)  shows  the  same  coher¬ 
ence  spectrum  for  the  Yj  stack  of  487  stations  outside  the  North  Sea  region.  There,  the  correlation 
between  pressure  and  sea  level  in  the  pole  tide  band  is  somewhat  below  the  95  per  cent  confidence 
level.  (There  is  a  large  coherence  at  a  period  of  about  460  days,  just  left  of  the  pole  tide  band.) 

Note  that  for  all  three  shallow  seas,  there  are  narrow  band  estimates  of  coherence  between  sea  level 
and  pressure  above  the  95  per  cent  confidence  level  at  several  frequencies  outside  the  pole  tide  band. 
While  it  is  not  unusual  for  the  oceans  to  be  stressed  by  winds  and  pressure  enough  to  show  enhanced 
coherence  at  a  particular  frequency,  what  is  coincidental  is  the  existence  of  a  normal  mode  of  the  earth 
(the  Chandler  wobble)  near  one  of  the  peaks  in  the  local  pressure  spectrum,  and  the  fact  that  the 
meteorological  forcing  in  this  case  is  concentrated  in  a  frequency  band  that  is  about  as  narrow  as  the 
pole  tide  peak. 

In  order  to  test  for  the  existence  of  a  high  Q  periodic  signal  in  the  atmosphere,  a  statistical  F-test  was 
calculated  from  the  Thomson  spectrum  of  the  pressure  signal  for  these  three  regions.  In  all  three 
regions,  the  F-test  value  is  greater  than  the  90  per  cent  confidence  level,  suggesting,  without  any  clues 
as  to  the  cause,  that  there  is  a  periodic  signal  in  the  atmosphere  this  region  with  a  period  of  14  months. 
The  F-test  exceeds  the  90  per  cent  confidence  level  in  several  other  frequency  bands  as  well. 

We  are  not  suggesting  that  this  is  evidence  of  a  pole  tide  in  the  atmosphere.  In  fact,  the  14  month 
pressure  peak  is  not  a  unique  feature  of  the  pressure  spectrum  (sec  Figure  VI. 5  (c),  below,  where  the  14 
month  peak  for  pressure  docs  not  stand  out  among  the  other  peaks). 

To  test  for  the  existence  of  a  resonance  response  of  the  shallow  basins,  the  sea  level  data  for  each 
basin  are  fit  in  the  frequency  domain  to  resonant  pressure  and  polar  motion  responses.  For  example,  let 
y(co)  be  the  Fourier  transform  of  the  stacked  data.  Since  wind  data  are  not  included  here,  we  attempt 
to  account  for  the  wind-driven  response  by  assuming  it  is  proportional  to  the  response  to  pressure. 
Then,  we  parameterize  y(cu)  using: 


y(o>)=  A  p(o))  +  B  r(co)  + 


C  p  (to) 

u  -  w0 


+ 


D  r(m) 
0)-  0)o 


+  E(CD) 


(7) 


Here,  A  and  B  arc  the  factors  that  describe  the  static  response  to  pressure  and  winds  and  to  the  polar 
motion  forcing,  respectively  (p(u))  and  r(co)  are  the  Fourier  transforms  of  the  pressure  and  equilibrium 
pole  tide);  C  and  D  represent  the  resonant  strengths  for  the  meteorological  forcing  and  the  pole  tide; 
co,  is  the  proposed  resonance  frequency  in  the  basin;  and  e(co)  represents  un modeled  effects,  including 
noise.  If  the  observed  peak  tn  the  basin  nca  the  pole  tide  frequency  band  is  due  to  wind  forcing  that  is 
correlated  with  pressure,  then  C  and  D  wiii  be  small,  and  y(u>)  will  be  dominated  by  the  static  terms. 
If  the  observed  peaks  are  due  to  a  resonance,  then  C  and  D  will  not  be  small,  and  to,  should,  presum¬ 
ably,  be  close  to  the  Chandler  Wobble  frequency. 

First,  we  find  the  values  of  A ,  B ,  C ,  D ,  and  cu0  that  minimize  e(co)  in  a  least  squares  sense.  The 
inversion  for  to,,  is  non-linear,  and  so  our  procedure  is  to  choose  a  great  many  values  for  co„  near  the 
pole  tide  frequency,  and  fit  for  A ,  B ,  C ,  and  D  for  each  of  these  co0 .  We  then  look  for  an  0)o  which 
leads  to  a  particularly  pronounced  decrease  in  the  residuals.  For  all  three  basins,  there  arc,  indeed,  pro¬ 
nounced  reductions  of  the  residuals  when  <a0  is  near  the  437  day  peak  in  the  pole  tide  band.  But,  wc 
judge  the  results  to  be  inconclusive,  due  to  the  high  correlation  between  the  two  resonant  terms  and 
B  r (m):  when  to,  is  close  to  the  pole  tide  period,  it  is  difficult  to  separate  the  resonant  term  from  the 
equilibrium  tide. 

Our  second  method  of  using  (7)  to  test  whether  the  data  require  an  oceanic  resonance,  is  to  simul¬ 
taneously  fit  and  remove  from  the  PSMSL  stacks  the  non-resonant  pressure  and  equilibrium  pole  tide 
(the  terms  proportional  to  A  and  B  in  cq.  7).  We  band-pass  the  data  to  extract  the  pole  tide  band,  prior 
to  the  fit.  Results  from  the  fits  suggest  that  the  pole  tide  in  these  shallow  seas  is  much  closer  to  equili¬ 
brium  than  is  suggested  by  simply  comparing  sea  level  with  the  equilibrium  tidal  amplitude.  For  exam¬ 
ple,  the  results  of  the  North  Sea  fit  imply  that  the  pole  tide  has  an  amplitude  of  1.2  ±  1.2  limes  the 
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equilibrium  amplitude,  and  that  the  response  to  meteorology  is  2.5  ±  1.3  times  the  inverted  barometer 
pressure.  For  the  Baltic  Sea  the  fit  to  the  equilibrium  pole  tide  is  1.9  ±2.1  fit  to  pressure  is  5.4  ±  1  8. 
For  the  Gulf  of  Bothnia  the  fitted  parameters  are  1.7  ±  1.0  and  4.5  ±  0.4.  These  large  formal  errors  on 
the  fits  arc  the  direct  result  of  large  variances  of  small  numbers  of  data  points.  The  pole  tide  frequency 
band  encompasses  only  5  bin  frequencies,  hence  the  two-way  fit  has  only  3  degrees  of  freedom.  Note 
that  none  of  the  pole  tide  parameters  arc  significantly  different  from  unity.  For  any  of  these  three 
basins,  if  the  fitted  pressure  and  pole  tide  are  removed  from  the  PSMSL  stack,  the  spectrum  of  the  resi¬ 
duals  (not  shown)  exhibits  none  of  the  characteristics  of  the  equilibrium  pole  tide  spectrum. 

To  help  assess  how  strongly  the  data  require  a  pole  tide  amplitude  greater  than  the  equilibrium 
amplitude,  we  subtract  the  equilibrium  pole  tide  from  the  data  (parameter  B  in  (7)  is  constrained  to  be 
unity),  and  then  fit  and  remove  the  pressure  (parameter  A,  alone)  from  the  PSMSL  stack.  If  the  spec¬ 
trum  of  the  residuals  show  the  characteristic  double  peak  of  the  pole  tide  spectrum,  then  there  is  reason 
to  suspect  an  amplified  response.  If  not,  then  the  fact  that  the  amplitude  of  B  exceeded  unity  in  the  pre¬ 
vious  fit  may  be  attributed  to  noise  leaking  into  that  fit. 

Figures  VI.5  (a)  and  (b)  show  amplitude  spectra  of  yj  stacks  for  the  Baltic  Sea  and  the  Gulf  of 
Bothnia.  The  sold  line  in  each  Figure  Vl.is  the  PSMSL  data.  The  spectra  of  the  PSMSL  stack  with 
pressure  fitted  and  removed  arc  the  short  dashed  curves  in  Figures  VI.5  (a)  and  (b).  The  long-dashed 
curves  are  the  spectra  of  the  residuals  after  subtracting  the  equilibrium  pole  tide  (short  dashed  curve) 
and  fitting  and  removing  pressure  from  the  PSMSL  stack.  The  spectra  show  that  the  removal  of  non- 
rcsonant  pressure  removes  the  majority  of  the  enhanced  power  at  the  pole  tide  frequency.  The  spectra 
of  the  residuals  after  subtracting  the  equilibrium  pole  tide  and  then  fitting  and  removing  pressure  have  a 
double  peak  but  they  do  not  follow  the  characteristic  peaks  of  the  equilibrium  pole  tide  closely.  The 
North  Sea  results  (not  shown)  arc  similar. 

The  problems  brought  about  by  the  lack  of  wind  data  in  this  analysis  are  illustrated  in  Figure  VI.5  (c) 
which  shows  power  spectra  of  Y\  stacks  of  pressure  (the  dashed  line)  and  PSMSL  data  (the  solid  line) 
over  the  Baltic  Sea  for  a  wide  frequency  band.  The  peak  at  the  pole  tide  frequency  in  the  pressure  spec¬ 
trum  is  evident,  but  it  is  not  appreciably  larger  than  other  pressure  peaks.  On  the  other  hand,  the  peak 
in  the  pole  tide  band  for  the  PSMSL  data  docs  stand  above  other  peaks  in  this  red  spectrum.  We  can 
also  infer  there  is  something  unusual  happening  in  the  pole  tide  band  from  the  least  square  fit 
coefficients  given  above.  For  example,  we  have  seen  that  the  fit  to  pressure  for  the  Baltic  Sea  gives  5.4 
±1.8  across  the  pole  tide  band  alone,  but  that  it  gives  1.2  ±  .1  if  the  data  are  not  band-passed  prior  to 
the  fit. 

Why  is  there  such  a  large  peak  in  sea  level  but  not  in  pressure  at  the  pole  tide  frequency?  There  arc 
several  possible  explanations.  It  could  be  that  there  is  a  broad  oceanic  resonance  across  this  band  that 
amplifies  the  response  to  the  meteorology.  Perhaps  the  wind  forcing,  but  not  the  pressure,  is  especially 
large  in  this  frequency  band.  Or,  maybe  the  wind  is  more  highly  correlated  with  pressure  in  this  band, 
so  that  the  wind-driven  and  pressure-driven  responses  tend  to  add  constructively. 

We  conclude  that,  overall,  our  results  arc  most  consistent  with  the  hypothesis  that  meteorological 
forcing  near  the  pole  fide  frequency  is  responsible  for  most  of  the  spectral  peak  in  sea  level  at  434-437 
days.  We  emphasize  that  the  PSMSL  data  that  lie  in  the  North  sea  region  are  limited  to  the  basin  per¬ 
imeters  and  have  small  numbers  of  stations,  that  we  did  not  examine  wind  data,  and  that  the  frequency 
band  of  interest  is  narrow.  We  are,  thus,  not  willing  to  conclusively  rule  out  resonance  as  a  mechanism 
for  the  enhanced  spectral  peak  in  these  basins. 


The  18.6  year  tide. 

Yearly  tide  gauge  data  arc  stacked  against  the  equilibrium  18.6  year  tide  for  all  260  stations  ouLside 
the  Baltic  Sea  that  reported  at  least  19  years  of  data.  The  Baltic  Sea  is  the  only  geographic  area  where 
long  period  noise  (there  is  a  spectral  peak  at  22-23  years  in  the  Baltic  Sea)  so  seriously  masks  the  18.6 
year  signal  as  to  affect  the  estimated  phase  of  a  global  stack.  We  identify  the  Baltic  Sea  as  a  particu¬ 
larly  noisy  region  as  follows.  First,  annual  and  18.6  year  amplitudes  are  fitted  and  removed  from  the 
data.  Then,  the  standard  deviation  for  the  stacked  data  set  is  calculated,  and  we  obtain  a=39.5  mm. 
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The  majority  of  stations  with  outliers  greater  than  3  (7  are  located  in  the  Baltic  Sea.  The  a  for  the  all 
stations  outside  the  Baltic  Sea  is  34.8  mm.  Data  outliers  are  not  removed  from  any  of  the  260  stations 
outside  the  Baltic  Sea,  as  they  do  not  affect  the  results  for  either  the  pole  tide  or  the  18.6  year  tide 
enough  to  justify  subjectively  altering  the  data  by  their  removal. 

Figure  VI.6  shows  PSMSL  and  synthetic  stacks  against  the  equilibrium  £(0,X)  in  the  time  domain.  A 
variation  of  approximately  16-20  years  may  be  seen  in  the  PSMSL  slack.  Figure  VI.7  (a)  shows  the 
amplitude  spectra  of  ungridded  equilibrium  stacks.  An  18.6  year  peak  stands  out  clearly  above  the 
noise  in  the  PSMSL  data,  and  has  an  amplitude  1 .2  times  the  synthetic  amplitude.  The  PSMSL  results 
also  show  pronounced  peaks  at  13.3  and  11  year  periods.  The  spectrum  of  the  difference  between  the 
the  synthetic  and  PSMSL  data  shows  that  most  of  the  discrepancy  between  the  PSMSL  and  synthetic 
data  at  the  18.6  year  period  may  be  attributed  to  noise,  as  the  amplitude  at  the  18.6  year  period  is  well 
below  the  noise  continuum.  The  results  suggest  that  the  PSMSL  and  synthetic  data  agree  well  in  both 
amplitude  and  phase  at  the  18.6  year  period.  Amplitude  spectra  for  a  stack  against  the  equilibrium 
C(0,X)  on  a  15°  by  30°  grid  is  shown  in  Figure  VI.7  (b).  The  use  of  a  coarse  grid  improves  the  agree¬ 
ment  between  the  PSMSL  and  synthetic  amplitude  at  18.6  years. 

A  statistical  F-test  of  the  multi-taper  spectrum  of  the  unweighted,  ungridded  data  is  shown  in  Figure 
VI.8.  (The  multi-taper  spectrum,  itself,  is  not  shown  here,  it  gives  results  similar  to  those  shown  in 
Figure  VI.7.)  The  F-test  shows  that  wc  can  be  95%  percent  confident  that  the  18.6  year  signal  of  the 
PSMSL  data  is  sinusoidal.  Note  that  the  13.3  and  11  year  peaks  are  not  as  prominent  in  the  F-test 
results. 

The  phase  of  the  difference  between  the  time  series  is  found  by  fitting  an  18.6  year  sinusoid  to  both 
the  PSMSL  and  the  synthetic  signals  and  comparing  their  phases.  The  uncertainty  for  the  amplitude  is 
estimated  as  described  above  for  the  pole  tide  results,  and  is  calculated  here  using  an  estimated  noise 
level  of  5  mm  for  both  the  gridded  and  ungridded  stacks.  The  results  are: 

Ungridded  PSMSL  is  1.18  ±  .22  times  the  synthetic; 

PSMSL  leads  by  1°±30° 

Gridded  PSMSL  is  1 .07  ±  .22  times  the  synthetic 
PSMSL  lags  by  10°  ±  30° 

We  combine  these  results  to  estimate  an  amplitude  for  the  18.6  year  tide  of  1.13  ±  .22  times  the 
equilibrium  tide,  and  the  phase  of  the  tide  to  be  equal  to  the  equilibrium  phase  within  a  30°  uncertainty. 
A  30°  phase  difference  corresponds  to  a  time  lead  or  lag  of  19  months. 

The  agreement  between  the  PSMSL  and  synthetic  18.6  year  spectral  peaks  is  not  overly  sensitive  to 
the  removal  of  blocks  of  stations  (with  the  exception  of  the  Baltic  Sea  stations,  which  have  already 
been  removed)  or  to  the  use  of  lime  spans  different  from  1900-1979.  For  example,  when  all  stations  in 
other  inland  seas  are  removed  from  the  data  set,  the  amplitude  at  18.6  years  is  found  to  be  affected  by 
an  amount  that  is  within  the  range  of  the  quoted  uncertainties,  although  the  frequency  content  of  the 
noise  does  change. 


The  global  rise  in  sea  level 


An  observed  secular  sea  level  change  at  an  individual  station  is  not,  by  itself,  evidence  of  a  global 
rise  in  sea  level.  There  could  also  be  secular  variations  due  to  post-glacial  rebound,  local  tectonic 
motion,  or  a  shift  in  the  wind-driven  oceanic  circulation  pattern.  The  effects  of  these  additional  secular 
changes  should  be  reduced  in  averages  of  global  data.  It  is  difficult,  though,  to  adequately  remove  the 
effects  of  post-glacial  rebound  by  averaging  alone  [Peltier,  1986).  Apparent  changes  in  sea  level  at 
individual  stations  due  to  post-glacial  rebound  can  be  as  large  as  8.5  mm/yr,  as  is  the  case  in  the  Gulf 
of  Bothnia.  And,  a  disproportionately  large  percentage  of  tide  gauges  arc  in  the  northern  hemisphere, 
close  to  the  centers  of  rebound.  Large  numbers  of  stations  also  lie  in  tectonically  active  areas,  and  no 
reliable  model  exists  that  allows  us  to  remove  the  crustal  motion  from  all  these  stations. 

Wc  digitized  the  rebound  results  from  Peltier’s  [1986]  post  glacial  rebound  models  for  North  Amer¬ 
ica  and  northern  Europe,  and  from  Wagner  and  McAdoo  [1986]  for  the  remainder  of  the  globe.  In 
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order  10  assess  global  sea  level  changes  most  accurately,  we  simultaneously  fit  the  entire  data  set  to  a 
linear  trend  and  to  an  synthetic  data  set  consisting  of  a  set  of  trends  predicted  for  each  station  by  the 
combined  rebound  models  as  described  above. 

In  order  to  establish  a  uniform  benchmark  for  each  station,  a  line  is  fit  to  each  station  and  the  y- 
interccpt  of  this  line  is  subtracted  from  each  station  record  before  the  simultaneous  fit.  The  intercepts 
subtracted  from  each  record  in  this  way  are  not  the  true  intercepts  for  those  stations  having  less  than  80 
years  of  data.  For  those  subsets  of  stations  containing  many  short  records,  this  introduces  a  systematic 
error  into  the  final  trend.  For  example,  if  a  station  record  having  only  a  few  years  of  data  contains 
power  at  periods  greater  than  the  record  length,  the  truncated  periodic  signals  are  correlated  with  the 
secular  variability,  and  any  true  secular  trend  over  several  decades  could  be  masked  or  even  reversed. 

On  the  other  hand,  if  we  restrict  the  data  set  to  only  those  stations  having  80  years  or  more  of  data, 
there  might  not  be  enough  stations  to  optimally  average  out  the  secular  variability  caused  by  sources 
other  than  the  global  sea  level  rise.  To  compromise,  we  first  construct  a  subset  of  the  data  that  includes 
all  stations  with  at  least  N  years  of  data,  where  N  is  greater  than  1  but  probably  less  than  80.  Wc 
choose  N  so  that  when  wc  simultaneously  fit  the  rebound  and  sea  level  rise  to  the  data,  the  post  glacial 
rebound  fit  parameter  is  close  to  1 . 

For  an  initial  data  set,  wc  choose  N=37.  This  leaves  us  with  120  stations,  and  the  simultaneous  fit 
results  in  a  rebound  coefficient  of  .9.  The  increase  in  sea  level  revealed  by  the  fit  is  1.2  ±  .1  mm/yr. 
The  uncertainty  is  the  rms  value  of  the  time  series  of  global,  yearly  averages  of  all  120  stations,  after 
the  rebound  is  removed.  A  post  glacial  coefficient  of  .9  is  a  good  indication  both  that  the  post  glacial 
models  are  giving  reasonable  results,  and  that  the  global  tide  gauge  data  are  capable  of  resolving  linear 
trends  on  the  order  of  millimeters  per  year. 

In  fact,  we  have  been  able  to  improve  the  fit  by  further  restricting  our  global  data  set.  First,  wc 
exclude  23  additional  stations  from  regions  of  the  globe  having  significant  tectonic  activity  (all  stations 
on  the  west  coast  of  North  America  and  Japan).  When  the  97  remaining  stations  in  this  reduced  data 
set  arc  simultaneously  fit  to  the  rebound  data  and  the  global  sea  level  rise,  the  rebound  coefficient  is  .94 
and  the  rise  in  sea  level  is  1.6  ±  .12  mm/ycar, 

Wc  further  exclude  all  stations  south  of  30°  north  latitude.  The  post  glacial  rebound  in  this  area  is 
small  and  is  reasonably  the  same  everywhere  (and  so  it  is  not  easily  separable  from  a  global  rise  in  sea 
level),  and  is  strongly  dependent  on  assumptions  about  the  Pleistocene  de-glaciation  of  Antarctica.  This 
final  data  set  consists  of  84  stations,  all  with  at  least  37  years  of  data,  situated  north  of  30°  N  latitude 
and  away  from  Japan  and  the  west  coast  of  North  America.  The  post  glacial  rebound  data  for  all  these 
stations  arc  estimated  from  the  results  of  Peltier’s  [1986]  model.  For  these  stations,  the  fit  to  the 
rebound  model  is  especially  good  at  .994  and  the  fit  to  the  linear  trend  is  1.75  ±  .13  mm/yr  over  the  last 
80  years.  This  trend  is  shown  in  Figure  VI.9. 

From  these  three  results,  wc  estimate  a  global  sea  level  rise  during  the  last  8  decades  of  between  1.1 
and  1.9  mm/yr,  with  a  preferred  value  near  1.75  mm/yr.  And,  because  the  post  glacial  fit  parameter  is 
so  close  to  1,  we  tentatively  conclude  that  the  effects  of  post-glacial  rebound  are  well  described  by 
current  models,  and  that  the  tide  gauge  data  arc  capable  of  resolving  global  linear  trends  at  the  millime¬ 
ter  per  year  level. 


CONCLUSIONS 

By  stacking  global  tide  gauge  data,  wc  can  significantly  improve  the  signal-lo-noise  ratio  for  long 
period  tides,  compared  with  results  obtained  from  single  station  records.  We  find  that  the  18.6  year  tide 
and  the  pole  tide  are  clearly  evident  above  the  noise  in  the  spectra  for  the  equilibrium  stacks.  The 
peaks  are  not  evident  in  stacks  against  other  pure  spherical  harmonics.  The  observed  amplitude  and 
phase  for  the  18.6  year  tide  are  consistent  with  the  assumption  of  global  equilibrium.  The  amplitude 
and  phase  of  the  pole  tide  arc  close  to  the  equilibrium  amplitude  and  phase  for  all  regions  outside  the 
North  Sea,  Baltic  Sea  and  Gulf  of  Bothnia.  In  these  basins,  the  enhanced  tide  appears  to  be  more 


consistent  with  narrow  band  meteorological  forcing  centered  about  the  pole  tide  frequency  than  to  a 
resonant  response  of  the  oceanic  basin,  although  the  results  are  not  entirely  conclusive.  We  find  that  in 
all  regions  except  the  Baltic  Sea  and  Gulf  of  Bothnia,  the  oceanic  response  to  atmospheric  pressure  are 
consistent  with  inverted  barometer  hypothesis  for  periods  longer  than  two  months.  Our  conclusions  are 
based  on  the  assumption  that  by  using  many  stations,  we  were  able  to  break  the  correlation  between 
winds  and  pressure.  This  cannot  be  confirmed  without  including  wind  data  in  our  analysis. 

Finally,  by  simultaneously  fiuing  the  global  the  global  sea  level  data  to  a  linear  trend  and  data 
prepared  from  a  combined  post-glacial  rebound  model,  we  conclude  that  the  global  rise  in  sea  level 
over  the  last  several  decades  was  between  1.1  and  1.9  mm/year,  with  a  preferred  value  near  1.75 
mm/yr.  Furthermore,  the  good  agreement  with  the  results  from  the  post  glacial  models  suggest  that 
those  results  are  reasonably  representative  of  the  true  uplift,  and  that  the  tide  gauge  data  are  capable  of 
resolving  global  changes  in  sea  level  at  the  millimeter  per  year  level. 
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CHAPTER  VI  FIGURE  CAPTIONS 


Figure  VI.  1.  The  coefficients,  A  (t),  of  the  least  squares  fit  of  sea  level  to  atmospheric  pressure  for  487 
stations  outside  the  North  and  Baltic  Seas  and  Gulf  of  Bothnia.  Seasonal  averages  have  been  removed 
from  both  the  pressure  and  sea  level  data.  The  response  is  consistent  with  the  inverted  barometer  assump¬ 
tion  for  time  lags  of  up  to  ±  240  months. 

Figure  VI.2  Amplitude  spectra  at  the  Chandler  wobble  frequency  band  for  Yj  stacks  of  487  stations. 
The  solid  curve  is  for  the  PSMSL  stack,  the  long  dashed  curve  is  for  the  synthetic  (equilibrium)  data  and 
the  short  dashed  curve  is  the  spectrum  obtained  when  the  synthetic  stack  is  subtracted  from  the  PSMSL 
stack.  The  effects  of  pressure  have  been  subtracted  from  the  PSMSL  data,  assuming  an  inverted  barometer 
response.  The  Chandler  wobble  is  clearly  evident  in  the  PSMSL  results,  and  its  amplitude  is  close  to 
equilibrium.  Stations  in  the  North  and  Baltic  Seas  and  in  the  Gulf  of  Bothnia  are  not  included  in  the  stacks. 

(a)  Ungridded.  (b)  10°  by  20°  grid. 

Figure  VI.3  Amplitude  spectra  of  Y  j  stacks  of  sea  level,  pressure,  and  the  equilibrium  pole  tide  (solid, 
short  dashed,  and  long  dashed  curves  respectively)  for:  (a)  the  North  Sea  (sea  level  is  3.6  times  the  equili¬ 
brium  pole  tide  at  the  437  day  period);  (b)  the  Baltic  Sea  (sea  level  is  6.2  times  equilibrium  at  the  437 
period);  (c)  the  Gulf  of  Bothnia  (sea  level  is  9  limes  equilibrium  at  the  437  day  period);  (d)  the  global  oce¬ 
ans  outside  the  North  Sea  region  (the  sea  level  and  equilibrium  spectra  are  in  good  agreement,  and  there  is 
no  prominent  pressure  peak  in  the  pole  tide  frequency  band).  This  figure  differs  from  Figure  VI.2  (b)  in 
that,  here,  seasonal  averages  have  been  removed  from  the  PSMSL  data  but  pressure  has  not.  The  pressure 
data  have  been  multiplied  by  10  mm/mbar  prior  to  constructing  the  stack.  Note  the  similarity  between  the 
sea  level  and  pressure  spectra  for  each  of  the  three  shallow  basins.  In  this  figure  and  in  Figure  VIA,  the 
peaks  in  the  pressure  and  sea  level  are  as  narrow  as  the  pole  tide  band  itself. 

Figure  VIA  Coherence  spectra  for  sea  level  versus  pressure  for  Yj  stacks  of:  (a)  26  North  Sea  stations; 

(b)  42  Baltic  Sea  stations;  (c)  19  stations  in  Gulf  of  Bothnia;  (d)  487  stations  for  global  oceans  outside  the 
North  Sea  region.  In  all  four  parts  of  this  figure,  seasonal  averages  have  been  removed  from  sea  level  and 
pressure  prior  to  constructing  the  coherence.  There  are  several  frequency  bands  where  the  coherence 
exceeds  .43  (the  95  per  cent  confidence  level),  including  the  pole  tide  frequency  band  for  the  three  shallow 
basins. 

Figure  VI.5  Amplitude  spectra  for  Y 2  stacks  from  stations  in  the  North  Sea  region,  (a)  Stacks  for  42  sta¬ 
tions  in  the  Baltic  Sea.  The  solid  curve  is  PSMSL.  The  pressure  is  fitted  and  removed  from  the  PSMSL 
data,  giving  the  short  dashed  curve.  The  equilibrium  tide  is  removed  and  the  the  pressure  is  fitted  and 
removed  to  give  the  long-dashed  curve.  Note  that  the  long-dashed  curve  does  not  follow  the  features  of  the 
equilibrium  pole  tide,  which  is  the  dash-dot  curve  near  the  bottom,  (b)  Same  as  (a),  but  for  the  19  stations 
in  the  Gulf  of  Bothnia,  (c)  Stacks  for  the  Baltic  Sea,  shown  over  a  large  frequency  range.  The  solid  curve 
is  the  PSMSL  data,  and  the  short  dashed  curve  is  the  pressure.  Both  spectra  show  peaks  in  the  pole  tide 
band,  but  the  PSMSL  peak  stands  above  all  other  spectral  peaks,  the  pressure  peak  does  not 

Figure  VI.6  An  ungridded  Y 2  stack  of  260  stations  as  a  function  of  time.  The  16-20  year  variation  in  the 
PSMSL  stack,  solid  curve,  is  discernible,  as  is  increased  noise  early  in  the  time  span.  The  dashed  curve  is 
the  stack  of  synthetic  data. 

Figure  VI.7  Amplitude  spectra  for  Y§  slacks  of  260  stations  outside  the  Baltic  sea.  The  solid  curve  is  the 
spectrum  for  the  PSMSL  data,  the  long  dashed  curve  is  for  the  synthetic  stack,  and  the  short  dashed  curve 
is  a  spectrum  of  the  difference  between  the  two  stacks  (a)  Ungridded.  (b)  15°  by  30°  grid. 
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Figure  VI.8  F-test  derived  from  Multi-taper  spectra  for  an  ungridded  Y§  stack  of  the  PSMSL  signal.  The 
F-tcst  predicts  with  a  95  percent  certainty  that  the  18.6  year  peak  is  sinusoidal. 

Figure  VI.9  Linear  trend  for  sea  level  for  84  stations  having  greater  than  37  years  of  data  that  lie  north  of 
30°  north  latitude  and  are  not  in  tectonically  active  areas  (Japan  and  the  west  coast  of  North  America  are 
excluded).  Vertical  motion  at  the  surface  is  simultaneously  fit  to  these  data  using  a  model  of  post  glacial 
rebound  from  Peltier  (1986).  The  fit  to  this  model  is  excellent  for  these  stations  at  .994,  and  the  fitted  rate 
of  global  increase  in  sea  level  is  1.75  ±  .13  mm/yr. 
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